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Abstract 

In this paper, on the basis of the Onsager-Wilson theory of strong binary electrolyte solutions we com- 
pletely work out the solutions of the governing equations (Onsager-Fuoss equations and Poisson equations) 
for nonequilibrium pair correlation functions and ionic potentials and the solutions for the Stokes equation 
for the velocity and pressure in the case of strong binary electrolyte solutions under the influence of an ex- 
ternal electric field of arbitrary strength. The solutions are calculated in the configuration space as functions 
of coordinates and reduced field strength. Thus the axial and transversal components of the velocity and 
the accompanying nonequilibrium pressure are explicitly obtained. Computation of velocity profiles makes 
it possible to visualize the movement and distortion of ion atmosphere under the influence of an external 
electric field. In particular, it facilitates tracking the movement of the center (x c , 0) of the ion atmoshphere 
along the x axis, as the field strength increases. Thus it is possible to imagine a spherical ion atmosphere 
with its center displaced to (x c ,0) from the origin. On the basis of this picture we are able to formulate a 
computation-based procedure to unambinguously select the values of x and r in the electrophoretic factor for 
£ > and thereby calculate the ionic conductance. This procedure facilitates to overcome the mathematical 
divergence difficulty inherent to the method used by Wilson in his unpublished dissertation on the ionic 
conductance theory (namely, the Onsager-Wilson theory) for strong binary electrolytes. We thereby define 
divergence-free electrophoretic and relaxation time factors which would enable us to calculate equivalent 
conductance of strong binary electrolytes subjected to an external electric field in excellent agreement with 
experiment. We also investigate the nature of approximations that yield Wilson's result from the exact 
divergence-free electrophoretic and relaxation time coefficients. In the sequels, the results obtained in this 
work are applied to study ionic conductivity and nonequilibrium pressure effects in electrolyte solutions. 
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I. INTRODUCTION 



Linear and nonlinear transport processes and nonequilibrium phenomena in dilute non-ionic 
(neutral) fluids have been known adequately treatable by means of singlet distribution functions 
obeying, for example, the Boltzmann equations and related kinetic equations[l| for singlet dis- 
tribution functions. Relying on the experience gained from the theories of neutral dilute fluids, 
theories of ionized gasesj^], plasmas, and charge carriers in semiconductors [3-0] often rely on singlet 
distribution functions obeying Boltzmann-like kinetic equations and their suitable modifications. 
However, since ions in ionized fluids interact through long-ranged Coulombic interactions, even if 
the ionized species are dilute in concentration, their spatial correlations are significant, lingering 
on to manifest their effects even in the infinitely dilute regime of concentration as the thermody- 
namic properties (e.g., activity coefficients) of ionic solutions demonstrate. Therefore it would be 
very important to find a way, and learn, to incorporate long-range correlations into the theory 
of nonequilibrium phenomena and transport processes in ionized fluids and therein lies the sig- 
nificance of the limiting theory of conductivity in ionized liquids in the external field of arbitrary 
strength described in this work. 

Interestingly, in the subject fields of nonlinear phenomena in ionic liquids, the Wien effect^] 
was one of the earliest experimental examples that exhibited a marked nonlinear deviation from 
the Coulombic law of conduction and, as such, it attracted considerable attention theoretically 
and experimentally. Being a nonlinear effect in ionic conductance which shows a strongly nonlin- 
ear, non-Coulombic field-dependence of ionic conductance, the phenomenon was studied actively 
in physical chemistry until several decades ago to understand ionic solutions and their physical 
properties ?[ sj]- Recently, there appears to be a revival of experimental studies on Wien effect and 
ated aspects in ionic conductance of ionic liquids in the presence of high external electric fields 



re 



11]. There are other many fascinating aspects of physica 
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properties of ionic liquids recently being 



13f ] . although they are mostly in the field 



studied actively and reported in the recent literature 
of equilibrium phenomena. In the present series of work, we are interested in nonlinear transport 
processes and, in particular, learning about the theories of the Wien effect on ionic conductance 
in electrolyte solutions in order to gain insights and theoretical approaches to treat the currently 
studied properties of ionic fluids. As a first step to this aim, we will study strong binary (symmet- 
ric) electrolytes because of the relative simplicity of the subject matter. More complicated systems 
of asymmetric electrolytes, in which the charges in a molecule are asymmetric, will be treated in 
the sequels [ul. to this work in preparation. 
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The ideas jjj3j] of physical mechanisms underlying the Wien effect, which might also encompass 
nonlinear phenomena in general in ionic fluids, proceed as follows. It is founded on the idea of ion 
atmosphere in Debye's theory fl?]] of electrolyte solutions. According to his theory, ion atmosphere 
is formed around ions in the solution, which is spherically symmetric if the ions are spherical 
and the system is in equilibrium. When the external electric field is applied to the ionic fluid, 
the ions of opposite charges begin to move in directions opposite to each other. Thus the basic 
physical mechanisms involved in the ionic movements under the external field are believed to be 
due to a distortion of the spherically symmetric ion atmosphere into a non-spherical form and its 
subsequent tendency to relax to a spherically symmetric form. The former effect gives rise to the 
electrophoretic effect and the latter to the relaxation time effect. It should be emphasized here that 
the aforementioned effects are on the ionic atmosphere, but not on the ion of attention situated at 
the center of ion atmosphere. 

This idea can be translated into a qualitative mathematical form as given below: In experiments, 
we measure migration of ions and accompanying flow of medium. If the external electric field is 
denoted X, the force kj on ion j of charge ej is then given by 

k i = e i x (j = V- - (!) 

Since the ion of charge ej in the solution creates an ion atmosphere of charge — ej, which is 
distributed in the ion atmosphere to balance the charge in the solution, and this atmosphere is 
subjected to a force of — e,X. This force tends to move the ion atmosphere in the direction of force 
— e^X, while the central ion j of atmosphere is carried by force e^X in the medium in the direction 
opposite to the motion of ion atmosphere in order to balance the momentum. The velocity of 
the countercurrent generated thereby may be readily calculated if it is assumed that the entire 
countercharge — ej of the atmosphere is distributed in a spherical shell of radius k -1 , where 
is the Debye radius of ion atmosphere from the central ion, and that the motion of this sphere 
of radius k _1 surrounding the central charge is governed by the Stokes law 1814211] holding for the 



motion of a sphere in a viscous fluid. Thus, this velocity of the countercurrent is estimated to be 

Av, = J — = — °- , (2) 

67T7/0 67T770 

where Avj is the velocity of the shell of radius kT 1 and r\§ is the viscosity of the medium. We are 
thus led to the result that the medium in the interior of the shell travels with this velocity, and 
that the central ion migrates against a collective current of the medium in the shell. The deduction 
of this expression qualitatively elucidates the most important part of the effect of electrophoresis. 
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Clearly, this effect has to do with hydro dynamic motion of the solvent around the center ion 
enclosed by the ion atmosphere of radius ft -1 that moves against the former. One may therefore 
quantify this qualitative description by means of a hydrodynamic method using the Navier-Stokes 



equation [19H2l|], but the Navier-Stokes equation requires a local body-force — local mean external 
force — as an input. This local body-force cannot be obtained through a purely phenomenological 
consideration, but, for example, must be calculated by means of statistical mechanics combined 
with classical electrodynamics. Before proceeding to the remaining effect, it is important to point 
out that Eq. ([2]) gives the velocity of a physical object of radius n~ l (i.e., the radius of ion 
atmosphere) in the direction of X. 

The second effect, that is, the relaxation time effect, is seen as follows: If the central ion 
possessed no atmosphere, it would simply migrate with a velocity where Q is the friction 

constant, but owing to its ion atmosphere, the ion is subjected to a net force, kj — Akj, where 
Akj is the force arising from the dissymmetry of the ion atmosphere created by the movement of 
the ions in the external field, and hence it will move, relative to its environment, with a velocity 
of a magnitude, (kj — Akj) /Q. This Akj is due to the effect arising from the relaxation of the 
asymmetric ion atmosphere. 



Consequently, the net velocity Vj of ion j is given by 

(e,X — Ak,-) e,-XK , N 

(3) 



Cj 6vr?7o 

Here Akj^" 1 represents the relaxation time effect on relaxation to a spherically symmetric form 
of the distorted ion atmosphere, and the last term the electrophoretic effect. 

The aforementioned two effects making up the velocity given in Eq. ([3]) are believed to underlie 
in charge conduction in electrolytic solutions. In fact, the mobility of ions induced by an external 
electric field can be calculated on the basis of the aforementioned two effects, for example, by using 
Eq. p. 

As we can see from this heuristic discussion, the aforementioned two effects require the velocity 
of the fluid (medium), which obviously obeys the hydrodynamic equations for the system subjected 
to an external electric field. Since such velocity solutions can be obtained from the Stokes equation, 
more generally, Navier-Stokes equation, we may apply the solutions thereof to calculate the charge 
conductance and the countercurrent of the medium to learn the mode of charge conductance in 
electrolyte solutions subjected to an external field. The hydrodynamic equations, however, contain 
external body-forces, which in the present case are the external electric field. The external electric 
field or body-force is generally local and depends on the local distribution of charges. The local 



charge distributions require molecular distributions in the system and a statistical mechanical 
theory for them — a molecular theory. 
To answer this question, Onsager 



22j with Fuoss formulated a formal framework of theory in 
which a Fokker-Planck-type differential equations for nonequilibrium pair distribution functions 
are derived on the assumption of a Brownian motion model for ions in a continuous medium of 
dielectric constant D and viscosity r/o- We will refer to these differential equations for pair correla- 
tion functions as the Onsager-Fuoss (OF) equations henceforth. They are coupled to the Poisson 
equations [2J] of classical electrodynamics for the ionic potentials. These two coupled systems of 
differential equations will be referred to as the governing equations in the present work. The 
governing equations were applied to study the ionic conductance of binary strong electrolytes in 
an external electric field by Wilson in his dissertation 25] . This theory will be referred to as the 



Onsager-Wilson (OW) theory. Wilson solved the governing equations and obtained analytic for- 
mulas for the electrophoretic and relaxation time coefficients and the equivalent ionic conductance 
qualitatively displaying the Wien effect in the regime of strong electric fields. Unfortunately, his 
dissertation has never been published in public domain, but only important results, such as the elec- 
trophoretic and relaxation time coefficients, had been excerpted in the well-known monograph^ 
by Harned and Owen on electro-physical chemistry. Tantalized by the possibility of the utility of 
the theory for recent experimental results for ionic fluids and charge carrier mobilities in semicon- 
ductors referred to earlier, we have thoroughly examined the OW theory to learn the details of it. 
Surprisingly, we have discovered that the velocity solution of the Stokes (hydrodynamic) equation 
in the OW theory can give rise to a divergent result rendering into question the electrophoretic 
coefficient calculated by Wilson's procedure described in his dissertation^25j. We believe that the 
basic framework of governing equations — the OF equations and Poisson equations — should be cor- 
rect, but the way the solutions are evaluated by him may be called into question. Therefore, it is 
our principal aim of this work to analyze the solutions of the governing equations in the case of 
binary strong electrolytes in an external electric field and obtain physically reasonable and thus 
acceptable theoretical results that can be made use of to study experimental data on conductivity 
and other transport phenomena in the high field regime. 

This paper is organized as follows. In Sec. II, we present the governing differential equations, 
which consist of the OF equations for the ionic pair distribution functions and the Poisson equations 
for the potentials of ionic interaction. We note that Kirkwood 26] also derived a similar equation for 



non-ionic liquids in his kinetic theory of liquids. One (BCE) of the present authors also derived 



the OF equations from the generalized Boltzmann equation. 



28, 



291 ] Since Wilson's dissertation 
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has not been published anywhere in a journal, the governing equations and their solutions are 
discussed to the extent that the present paper can be followed intelligibly. 

In Sec. Ill, the solutions of the governing equations — the pair distribution functions and poten- 
tials of ionic interaction — are presented in the case of a strong binary electrolyte solution subjected 
to an external field. These solutions are given in one-dimensional Fourier transforms in an axially 
symmetric coordinate system, namely, a cylindrical coordinate system whose axial coordinate is 
parallel to the applied external electric field. The Fourier transform is with respect to the axial 
coordinate. The distribution functions obtained are nonequilibrium pair distribution functions 
which describe the nonequilibrium ionic liquid structure, and the nonequilibrium ionic potentials 
of interaction in the external field. Since they should be of considerable interest to help us learn 
about the nonequilibrium ionic liquid properties we study the solutions of the governing equations 
in detail and obtain, especially, their spatial profiles, indicating how ions and their nonequilibrium 
part of the potentials are distributed in the external electric field. It should be noted that the 
distribution functions are the nonequilibrium corrections to the Boltzmann distribution function 
predicted by the Debye-Hiickel theory 17|] of electrolytes, and similarly for the potentials. 

In Sec. IV, we then discuss the solutions of the Stokes equation, which replaces the Navier- 
Stokes equation in the case of incompressible fluids that we assume the ionic solution of interest 
is. Solving the Stokes equation, we obtain the axial and transversal velocity components as well 
as the nonequilibrium pressure from the solutions of the Stokes equation. We present the solution 
procedure for the Stokes equation in detail, because, firstly, Wilson's thesis contains only the 
symmetric part of the solution, leaving out the antisymmetric part that turns out to be comparable 
to the former in magnitude and, secondly, we believe that the solution procedure of the Stokes 
equations, which combines statistical mechanics and hydrodynamics in a rather intriguing manner, 
appears to be very much worth learning, especially, if one is interested in nonequilibrium theories 
of ionic liquids in an external electric field. In this section we also discuss the connection with 
the electrophoretic and relaxation time coefficients originally obtained by Wilson, who evaluated 
them at the position of the center ion of ion atmosphere, namely, at the coordinate origin. This 
discussion would show that one of his integrals evaluated at the coordinate origin is divergent. 
Therefore we evaluate explicitly the solutions to explore a way to make the OW theory of ionic 
conductance unencumbered by such a divergence difficulty. 

In Sec. IV, we also compute numerically the spatial profiles of the axial velocity, and study 
them to guide us to avoid the divergence difficulty mentioned in connection with Wilson's re- 
sult and choose the optimum position coordinates at which to calculate the relaxation time and 
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electrophoretic coefficients. To this aim we have either evaluated analytically or reduced to one- 
dimensional quadratures, by means of contour integration methods, the Fourier transform integrals 
making up the solutions of the Stokes equations obtained earlier before computing their spatial 
profiles. The contour integration methods are described in Appendix A. Since they, however, do 
not cover the entire coordinate space owing to the condition imposed by Jordan's lemma 30(] on 



applicability of contour integration methods involving integrations along a circle of infinite radius, 
the integrals must be numerically computed outside the region where the aforementioned condi- 
tion is violated. The details of the condition are discussed in Sec. IV and also in Appendix A. 
These numerical studies reveal the manner in which the ions flow subject to the applied external 
electric field provide insight into the behavior of the velocity and valuable clues to formulate an 
empirical rule to select the position parameters (x, r) in the electrophoretic factor, so that a phys- 
ically sensible and non-divergent electrophoretic coefficient and the corresponding relaxation time 
coefficient can be defined and ionic conductance correctly predicted. This problem is addressed in 
the companion paper. Sec. V is for discussion and concluding remark. 

II. GOVERNING EQUATIONS 

Let Tj denote the position vector of ion j in a fixed coordinate system and Tji the relative 
coordinate of ion i from ion j: 

V = Tj{ — Fj Tj — Tj Fj — Tij = r, (4) 

and let fji (Tj,Tij) denote the concentration of ion i in the atmosphere of ion j located at position 
Tj — in other words, the distribution function to find ion i at distance from ion j located at Tj. 
At equilibrium it is given by the Boltzmann distribution function times the density of ion j. Let 
us also denote by Vji the velocity of ion i in the neighborhood of ion j. Therefore this velocity also 
depends on positions of ions i and j in the following manner: 

Vjj = Vji (Tj,Tij) , Vij = Vij (Ti,Tji) . (5) 

The equation of continuity for ion pair ( j, i) is then given by 

" i)J > ; *f' j] = V 3 ■ (vijfij) + V, • (v*/*) = , (6) 

where Vj = d/dfj. Hence, at a steady state dfji/dt = the steady-state equation of continuity is 
given by 

V ; ••:v, ; /, ; ! • V, -iv ; ,/ ; , ) 0. (7) 



Assuming that the ions, being randomly bombarded by molecules of the continuous medium (sol- 
vent) of dielectric constant D and viscosity 770, move randomly, namely, execute random Brownian 
motions, in the presence of an applied external field, the velocities Vjj and Vjj may be assumed 
given by the Brownian motion model 

Vji = V(fj) + Ui(Kji - k B TVi In fj i} (8) 
Vjj = V(r,) + Uj(Kij - k B TX7j In f ij} (9) 

where V(r&) is the velocity of solution at position {k = LOk is the inverse of the friction 
coefficient C/c of ion k, which is related to the diffusion coefficient D k of ion k of charge e k in the 
medium of viscosity rjo 

D k = k B Too k = ^L. (10) 

Here k B is the Boltzmann constant and T the absolute temperature; Kjj is the total force acting 
on ion i. We assume that forces on ions Kjj are linear with respect to charge numbers 

Kji = kj - eiViipj (r_j-,ry) , (11) 

so that the superposition principle of fields is preserved. Here kj is the applied external force on ion 
i. Under the assumptions for v™ and for Kj, stated earlier, the steady-state equation of continuity 



22f | satisfied by ion pair distribution functions 



((7J) becomes a coupled set of differential equations 
fji (rj,fij) of the ionic liquid: 

k B T (uji + ujj) V • Vfji (r) + (ajjCj - ^ej) X • Vfe (r) 

+ni nj [emJiV ■ (r) + ejujjV ■ V^i (-r)] = 0, (12) 

= 1,2, • • • ,s). 

We will call this set of differential equations the Onsager-Fuoss (OF) equations. Here for simplicity 
of notation we have omitted the first position variables in the distribution functions and potentials 
and typeset them as follows: fji(r) = fji (rj, r^) , etc. and ifij (r) = ipj(Yj,r) and tpi(-f) = 
tpi (Fj, — r). In fact, for Eq. (|12p the coordinate origin may be regarded as fixed on position of ion 
j. These are Fokker-Planck-type equations for fji (r) and tpj (r) and (— r). In Eq. (|12p . rij is 
density of ion i and X is the external (electric) field. The potentials ipk (k = appearing in this 



set of differential equations, Eq. (|12p . obey the Poisson equations of classical electrodynamics 



4-7T 

V-V^(r) = -— -^e^(r). (13) 



Dm 

i=l 



The two sets (|12p and (|13p are coupled to each other and will be henceforth referred to as the 
governing equations in this work. 



A. Boundary Conditions 

The two sets of equations (|12j) and (|13|) are subject to the boundary conditions stated below. 
1. No Flux Conditions 

The number of ions, leaving and entering the interior, f2, of a surface S should be balanced, 
because no ions are created or destroyed. Therefore fji(r) [vj{ (r) — (— r)] is sourceless in Q. 
This fact may be expressed as 

I dSfjiit) {e n ■ [ Vi< (r) - Vij (-r)]} = [ cMV ■ {/^(f) [v* (r) - v y (-f)]} = 0, (14) 

where e n is the vector normal to the surface S. This will be henceforth called no flux condition. 



2. Boundary Conditions on Potentials 

If the charge ej is within fi, we obtain 

lim / dVtpj (r) = e,-<$, (15) 
where (r) is the charge density. Therefore, the space charge within must be such that 

- ^- I dnV 2 ^j (r) = / dn Qj (r) , (16) 

or alternatively 

lim / dS ■ Vifjj (r) = -2-6, (17) 

for the boundary condition on the ionic potentials. Here 



8 



1 if ion j is located at r = 
otherwise 

The boundary condition (|1T|) corresponds to the fact that the charge ej at the origin (i.e., at the 
center of the ion atmosphere) must balance the net charge of the rest of the ion atmosphere. 

B. Symmetric and Antisymmetric Parts of Governing Equations 

Since it is convenient to work with dimensionless variables, we first reduce position variable r 
with respect to the Debye parameter k 

r = «x, (18) 



where the Debye parameter is defined by 



4vrr r 



ro = E 



n k e\. 



(19) 



k=l 



Dk B T 7 

The distribution functions and potentials change the sign of argument if particle indices j and i 
are interchanged. Therefore they are expected to consist of symmetric and antisymmetric compo- 
nents. Consequently, it is convenient to distinguish the symmetric and antisymmetric components 
/•£ (r) and tp- (r), etc. of distribution functions and potentials. They are defined in reduced forms 
as follows: 

1 



Dk B T Jji { 



[fji(r) + fji(~ r )} ~ n i n j. 



2 

J[^-(r)±Vi(-r)]. 



(20) 



(21) 



Since the distribution functions tend to riiUj as r = |r| tends to infinity, (r) vanishes as r — > oo. 
According to the definitions (|20f) and (|21fl . the following symmetry properties can be deduced for 
them as the ion positions are interchanged: 



fji (-*•) = fti w = 4 (- r ) - fji (- r ) = -fn w = -fii (- r ) > 



-r 3 (r) 



(22) 
(23) 



The differential equations of the symmetric and antisymmetric components of fji and ipj in Eqs. 
(fl~2l) and (fl"3j) then can be separated as follows: 



^ 2 fji (r) - (ujijfii + WjiMi) /J (r) - ^i.j I '.,/,) ( r ) - WjilHfu ( r ) + I'ji^-rf j; ( r ) = 0, 

V 2 // ; ir! y/;/ ?; !r) • //,/ ;j ir) 
V 2 /+(r)- [^(r)+^/+(r) 

V 2 /^ (r) - (wy^ + Uji/ij) (r) - u^^fr. (r) + w^/:: (r) + //',V.,./y, (r) 

V 2 /-. (r) 



0, 

o, 

0, 

o, 



(24) 



V 2 /^ (r) = 0, 



V V/ (r) 
V (r) 



V 2 V7 (r) 



v 2 vr to 



Mi/* W + to 
Mi/+(r) + Mi/+(r) 
/'//;,- to + /'.//,., (f) 
Mi/^ to + Mi/,^ to 



(25) 
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where V 2 now stands for Laplacian operator of reduced variable r, 

= 2~ 2> = , ' ( 26 ) 

u'aa = ^i 6 ! — u%e "i) (X = external electric field strength) , (27) 

n% nk B T ujj + uji K ' y ' 

with Vi denoting the stoichiometric coefficient of ion i in the (j, i) electrolyte and the charge 
number of ion k: = ezt (k = Henceforth the indices j and i refer to ions of the binary 
electrolyte (j, i), but also may dually refer to other ions belonging to species j or i. This nota- 
tional device prevents proliferation of subscripts distinguishing ionic particles. The 10 differential 
equations of Eqs. (|24p and (|25p will be referred to as the governing equations, the set (|24p as the 
Onsager-Fuoss (OF) equations, and the set ([25]) as the Poisson equations. The solutions of the 
governing equations provide the information on the nonequilibrium ionic liquid structure and ionic 
potentials of the electrolyte solutions subjected to an external electric field of arbitrary strength. 
A theory of transport processes in ionic solutions can be developed by making use of them. 

III. NONEQUILIBRIUM IONIC LIQUID STRUCTURE AND POTENTIALS OF BI- 
NARY ELECTROLYTES 

A. Complete Solutions of the Governing Equations 

n 

We now limit our study to strong binary electrolyte solutions as in the theory of Wilson [25] and 
Onsager. If the electrolyte is binary and strong, then \ej\ = \ei\ = ez with I, and 



1 , zeX 

Wilson in his unpublished PhD thesis [25f] obtained formal solutions of the governing equations in 
the forms of Fourier transforms under the assumption that 0J{ = ujj, which means ujji = uJij = \; 
that is, the diffusivities of the constituent ions are equal. (As it will turn out, the difference in the 
diffusivities has only a minor effect that can be ignored without much effect on the solutions.) And 
therewith he formulated a theory of ionic conductance of binary electrolytes under the influence of 
external electric field. However, Wilson's thesis unfortunately has not been published in a journal 
in public domain, nor have the nonequilibrium ionic liquid structures and accompanying poten- 
tials been explicitly evaluated and studied. In fact, neither were the velocity profiles completely 
calculated in the full configuration space since he limited the study to the velocity of the center 
ion of the ion atmosphere located at the coordinate origin in his calculation of the electrophoretic 
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effect. Moreover, the particular velocity formula made use of by Wilson had a divergence difficulty 
at the origin, but he argued it away on the ground that the divergent term would not contribute to 
the ionic conductance. We will show his argument was mathematically groundless and would not 
hold true. For these reasons, in this work we will first evaluate the velocity formulas explicitly by 
applying analytic methods or methods of contour integrations or numerical computation methods 
for wide ranges of position coordinates, and then will explore a way to overcome or get around the 
divergence difficulty. The results obtained thereby for the nonequilibrium ionic liquid structure 
and potentials as well as the velocity profiles would be principal contributions of the present work, 
which are not available in the literature on ionic liquids at present. In the subsequent paper 3l| . 



the solutions of the present paper will be applied to study the Wien effect on equivalent ionic 
conductance as a function of the applied field strength. 

Since Wilson's dissertation is not only not readily accessible as mentioned earlier, but also his 
solution procedure is difficult to follow, on the basis of our understanding of his solution procedure 
we will reconstruct the solutions for the governing equations (|24f) and (|25p . The solution procedure 
presented below is not exactly the same as his except in spirit, but most of the final results agree 
with his in the main. Under the assumptions on ojji = cj™ mentioned earlier, the governing 
equations (|24p and ([25]) are given as follows: 



V 2 f+ 



V 2 /^ W " i # W + ^ (r) 



1 



1 



V 2 4- (r) 
V 2 ^ M 



0, 

o, 
o, 
o, 
o, 

0, 



(28) 





(r) 


1 

~ ~2 


/* W 


+fii 


(r) 




(r) 


1 

~ ~2 




+$ 


(r) 




(r) 


1 

~ ~2 


/, w 


+/-• 


(r) 




(r) 


1 

~ ~2 






(r) 



(29) 



Owing to the fact that the solution of the Laplace equation is constant, the solutions of the 
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fifth and sixth equations of the set (|28|) are constant: 



fn W = C jt ffi (r) = C h 

but by the boundary conditions that they must vanish as r — > oo. Therefore, the constants Cj and 
Ci must be equal to zero. Hence 



4-(r) = 0, fu (r) = 0. 



(30) 



Consequently, the governing equations reduce to the following 8 differential equations: 



V 2 -^)/+(r)-^[/5(r) + /±(r) + (V,/-(r): .). 



V 2 " \) $ (r) ~ \f± (r) = 0, 
(v 2 - ^ (r) - \f± (r) = 0, 
V 2 -^/n( r ) + eV,M(r) = 0, 



2 / ^ 



(31) 
(32) 
(33) 
(34) 



V 2 ^ 



V 2 ^ 



v 2 ^- 



v 2 ^ 



|[/iW + /i(r) 



+ (r 



l r 

2 

2^?* ( r ) ' 
\fii (r) • 



(35) 
(36) 
(37) 
(38) 



These two sets, (|3Tj) - (|34|) and (j35]) - (j38|) . suggest that having obtained the solutions of the first 
set (|3T1) - (f3"4"j) we can look for the solutions of the second set (|35l) - (f38|) . inhomogeneous differential 
equations . We will follow this strategy by applying the method of Fourier transform. 

Since there exists an axial symmetry present in the system owing to the fact that a uniform 
external electric field is applied in a direction, we choose a cylindrical coordinate system whose 
axial coordinate axis is parallel to the external field direction. The cylindrical coordinates will be 
denoted (x, p, 6) where x is the axial coordinate, p the radial coordinate transversal to the axis x, 
and 9 the azimuthal angle; see Fig. 1. Then the distribution functions and potentials have axial 
symmetry around the x axis, and hence they are independent of angle 9. Now Fourier transforms 
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are taken with respect to the axial coordinate x: 
( 




7T 



da 



' cos (ax) f± (a, p) 



etc. 



^ Jo 



sin (ax) f j{ (a, p) 



(39) 



da 



etc. 



(40) 



cos (ax) i/jj (a, p) 
, sin (ax) ipj (a, p) 

Here a is a dimensionless wave number in units of k. When Fourier transformed in this manner, 
the governing equations (|3~T]) - (}3"4"]) and (f35]) - (f38j) become sets of coupled second-order ordinary 
differential equations with respect to the reduced radial coordinate p (perpendicular to the x axis) 
given below: 



l r 



^ ) fji (a, P) ~ 1 [fjj (a. P) + fu p) + of /,•< (a, p) = 



'" 2 
1 



) 2 

' 2 



2" 
1 



1 



/ji («. p) - a £/ji ( a > P) = °. 



(41) 
(42) 
(43) 
(44) 



Dffi (r) 
Dffi (r) 
£>?#7 (r) 
V% (r) 



§ (r) + f u (r) 
f+(r) + f±(r) 



1 

2 
1 

2 

2^ ( r ) ' 
'2~S W • 



Here symbol .D 2 is defined by the differential operator 

p pdp dp 



(45) 
(46) 
(47) 
(48) 

(49) 



Because the zeroth-order Bessel function Kq (ap) of second kind is an irregular solution of the 



differential equation 



32, 3] 



DlK Q (ap) 



Id d 

- p Tp p dp- a ) k o(«p) = o, 



(50) 



the coupled inhomogeneous differential equations (|4ip444p are expected to be solved by linear com- 
binations of zeroth-order Bessel functions but of different arguments XkP, where (k = 1,2,3,4) 
are characteristic values of the differential equation system. Unfortunately, two of the character- 
istic values turn out to be degenerate. Therefore it is not possible to apply the method of linear 
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algebra to solve the system in the conventional manner in which the solutions are expanded in 
characteristic vectors. This difficulty is overcome if Eqs. (|4ip - (|44|) are solved in the following 
manner. 

Operating [D 2 — |) on Eq. (jlTJ) and eliminating resulting (D 2 — i) ft, [D 2 — i) f£, and 
(D 2 p — ^) using Eqs. (j4"2j) - (l4"4"l) . we obtain the fourth-order differential equation 

(D 2 p - \ + ±Rj (d 2 - 1 -- ±Rj f+ (a, p) = 0, (51) 

where 



i?=Vl"4K). (52) 

This fourth-order differential equations can be solved by Bessel function Kq (Aip) and Kg (A 2 p), 
where Ai and A2 are two characteristic values 



A, M /I +a 2 + iji, \ 2 = \j l - + a 2 - l -R. (53) 



These are non-degenerate. Therefore the general solution for Eq. (|51l) may be written as a linear 
combination of the Bessel functions 

f± (a, p) = A l K (X lP ) + A 2 K (X 2 p) , (54) 

where A\ and A 2 are constant coefficients that must be determined by the boundary conditions, 
Eqs. (|14p and (|17p . or the equivalent conditions, for the symmetric and antisymmetric parts. Note 
that this solution satisfies the boundary condition as p — > 00 since the Bessel functions Kq {\kp) 
vanish at p = 00. Upon substituting this expansion into Eqs. (|4ip - (|44p we obtain 

-\ffj ^ ~ \f* ( a "°) + a ^Ji ^ = ~\ RA ^ ^P) + ( A ?p) , (55) 

D p ~ I) ft ( a ' P) = ( A ^) + \ A2K * ^ > ( 56 ) 



D 2 P ~ I) ft («,p) = \mK q (A 1P ) + ^A 2 if (A 2 p) , (57! 



D p " \) ft ( a ' ^) = a ^ K o (AlP) + a^2^o (A 2 p) • (58) 



This inhomogeneous set may be also solved by expansion. Let 



f± (a, p) = B 1 K (Aip) + # 2 ^ (MP) + B 3 K (A 3 p) , 

ft (a, p) = Cii^o (Aip) + C 2 if (A 2 p) + C 3 K (A 3 p) , (59) 
( a , p) = (Aip) + D 2 K (A 2 p) + D 3 K (A 3 p) , 



15 



where Bk, Cfc, and are expansion coefficients to be determined and A is the degenerate charac- 
teristic value to be determined self-consistently. Inserting these expansions into Eqs. (|55p - (|58|) we 
find relations between the coefficients and also the as- yet undetermined characteristic value A. We 
find 



A = y« 2 + ^ = A 3 (60) 

which is the degenerate third characteristic value of the governing OF equations for binary elec- 
trolytes. It is independent of the external field strength X or £ unlike Ai and A2. The relations 
between the coefficients are also obtained as follows: 

Bi = B 2 = - 2 -^A 2 , 

R R 

^1 = ^1, C 2 = ~A 2 , (61) 

Dl = I" 41 ' ° 2 = ~\ A ^ 

Thus the distribution functions f^, f£, are given as linear combinations of Bessel functions 
K (X kP ) (k = 1,2,3): 

f+ (a, p) = A X K (Aip) + A 2 K (X 2 p) , (62) 

S ^ = ^JT AlK ° (Al/a) " ^ A2 ^° (A2P) + BsK ° (Ap) ' (63) 

fu («, P) = ^i#o (Aip) - -^A 2 ^o (A 2 p) + C 3 ^ (Ap) , (64) 

ft (a, p) = i^iTo (Aip) - ^A 2 K (A 2 p) + D 3 ^ (Ap) . (65) 



The solutions of Poisson equations f|45|) — (|48[) can be similarly obtained as linear combinations of 
Bessel functions Ko^X^p): 

^ (a, p) = —AiKo (Aip) + ^A 2 K (A 2 p) - C 3 K (A 3 p) , (66) 

& (a, P) = -^Kq (Aip) + i^ 2 Jfo (A 2 p) - (4a££ 3 - C 3 ) J*T (A 3 p) , (67) 

(«, p) = - rK$r\ AiKo {Xlp) + Rn^R) AiKo {X2P) ~ BiKo (Asp) ' (68) 

$r («> p) = - R ?il R ) AiK * ^ + RTrhij A2Ko {X2P) " B3 ^ (A3P) • (69) 

The coefficients A\, A 2 , B3, and C 3 in these expansions are determined by imposing the boundary 
conditions (|14p and (|17fl . which for the symmetric and antisymmetric parts become 

s(^s + H^ + H* f )= ' (70) 

iT.(4 /i + "^) =0 <* = Ai) ' (71) 
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]jmp^+(p) = -6 (k = i,j), 
p^o dp 

]imp—^(p) = (k = i,j). 
p->o dp 

We note that the behavior of .Ko(Afcp) near p = has the property 



(72) 
(73) 



p— >o dp 



d 

lim p— lnp = — 1. 

p^o r dp 



(74) 



Imposing the boundary conditions, we obtain the linear algebraic relations of coefficients, which 
can be easily solved upon reducing them to independent linear equations. To save the space we 
simply present the final results only: 

(R + l) A (12-1) 



Bi 
Ci 
Di 



2R ' 
(R + 1) a£ 

R 2 
(R+l) 



2R 2 ' 
2R 2 ' 



A 2 

C 2 
D 2 



Bo 



2R ' 
(1 - R) a£ 



R 2 



2a£ 



2i? 2 ' 

2« 2 ' 



Da 



16 (a£)' 
i? 2 

4(aQ 2 
i? 2 



(75) 
(76) 
(77) 
(78) 



The Fourier components and in Eqs. (|85p and (|86p are finally given by 



# («. = [(« + !) (Aip) + (R - 1) i^o (A 2 p)] , 



2R 



/« («^) = + !) (Aip) + (fl- 1) (A 2 p) - 2K (Ap)] , 



'33 



2i? 2 



(-R + 1) (Aip) - (fl - 1) JKb (A 2 p) - 8 «) 2 ifo (Ap) 



f» («. = -^2 P + !) (Aip) + (1 — -R) ET (A 2 p) - 2K (A 3 p)] , 



2i? 2 



(i? + 1) K (Aip) + (1 - R) K (A 2 p) - 8 Kf (A 3 p) 



(«, P) = ^2 [#o (Aip) + #o (A 2 p) - 2K (A 3 p)] (fc = j, i) . 



(79) 
(80) 
(81) 
(82) 
(83) 
(84) 



We summarize the Fourier transforms of the solutions for the distribution functions and potentials 
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we set out to find: 

-2^ „ „ _ /-oo 



(r) - = ^ da — [(i? + 1) K (X lP ) + (R - 1) K (A 2 p)] 



'o 

2K 2 njriiejei^ f°° a sin (ax 



, da ^ — -x (85) 

[(R + 1) K (A lP ) + (1 - i?) JKb (A 2 p) - 2K (A 3 p)] , 



e k f°° , cos (ax) 
^k(r) = — da p2 x 



2D ./n i? 2 

(i? + 1) K (Aip) + (1 - i?) iT (A 2 p) - 8 (a£) 2 K (Agp) 

and similarly for fjj(x,p) and fu(x,p) to fji(x,p) in Eq. ([85]) . It should be noted that the 
variables in the integrals are reduced variables in the units of the Debye parameter k; see Eq. (|87l) 
below. 

The solutions presented in Eqs. (I85D and (I86p are the nonequilibrium parts of the pair distri- 
bution functions and potentials in the ionic liquid in the external field X (or £ in reduced form) 
at arbitrary strength. Therefore they represent the nonequilibrium ionic liquid structure and po- 
tentials when the ions are moving subjected to the external field at a steady-state condition. If 
the full potential is desired, ipf. (a, p) must be combined with the equilibrium Debye potential — the 
Yukawa-type potential. Therefore it would be of great interest to see how the nonequilibrium liq- 
uid structure and potentials vary with respect to spatial positions and the field strength. We will 
investigate these aspects (i.e., profiles) in the following. 



B. Evaluation of Nonequilibrium Ionic Liquid Structure and Potentials 



The Fourier integrals in Eqs. (I85h and (|86l) contain three parameters, position coordinates 
x and p and the reduce field strength £. Although looking complicated, they can be evaluated 
analytically in the region where the transversal (radial) coordinate p satisfies a certain condition 
with respect to the axial coordinate x, as will be stated more precisely later; see Eq. (190p below. 
In the rest of the (x, p) plane where the condition is not met, they can be evaluated numerically, 
provided that the singular behavior of the integrands is properly handled by using the method of 



principal values used for singular integrals 



34] 



For the evaluation of the integrals, it is convenient to scale further the variables as follows: 

t = V2a, r = p/V2, x = x/V2, u k = V2X k . (87) 
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Thus 

Ul = ^/i + t 2 + x /T3^2 5 U2 = ^l + ^.yr^W, to 3 = Vl + t 2 . 
It is also convenient to define 



001 = v i - y 2 + \/i + 2eV, zj 2 = V i - y 2 - \A + 2£V 

& 8 = >/r=v, a=v^3!. (89) 

Notice that = ^fc|t=iy (k = 1, • • • ) with the complex variable £ taken along the imaginary axis. 



As shown in Appendix A, if the condition 



30] 



x ReiOk(t) , x 

- + T — > (90) 
r Imi 



for x, r > in the complex plane of variable t, the integrals in Eqs. (J85J) and <\86h can be evaluated 
analytically or reduced to simple one-dimensional quadratures, if methods of contour integration are 
employed. As a matter of fact, the one-dimensional quadratures thus obtained can be analytically 
evaluated term by term in series if the series representation for Bessel functions Iq(z) is used. 
Condition (|90p means that the region in question is roughly within a conical domain surrounding 
the x axis. Outside this region the integrals must be computed numerically by applying methods of 
principal integrations for singular integrals 



34] . In this exterior region the integrals vanish uniformly 

88 a, ^ d , - — ^ „ . ,, « ^ ^ . 

the special position of x = r = 0, namely, the coordinate origin. Henceforth for notational brevity 
the reduced variable x will be simply typeset x without the caret^. 

The integrals in Eqs. (|85|) and (f86]) . reduced as described above, are evaluated by means of the 
contour integration methods described in Appendix A. They are given by the expressions 

* ( ±r ) = - 2 - u dye - ( 1 + TTTw J /o <^> 



/•x/^l+O e-^y ( 1 + X /1 + 2£V) 
^r^^^/oCaJsr)}, (91) 



o 
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$3 ( ±r ) = (^ r ) 



8ttV2D 



V^<J+F) e-*y (1 + yiT2£V 



dy 1 i oc2 2 7 ° ( Wir ) 

1 + 2£^2r 



f 1 P - x v v 2 



±\/2£ 



^ 2{1+?2) . ye'^/o far) /- 1 ye-^Io (w 3 r) 

dy—r—iTpr^ — 2 / 



(92) 



Here Iq facr) (k = 1,3) are the regular Bessel functions of zeroth order of second kind. In these 



expressions the range of position variables x and r must be such that -\/2 (1 + £, 2 )x > r for the inte- 
grals involving uJi, and x > r for the integrals involving TJ3. These conditions, related to Condition 
(|90p stemming from Jordan's lemma 30] on contour integrals, ensure the boundary conditions for 
the distribution functions and the potentials, which vanish as x and r tend to infinity. 

The results presented in Eqs. (|9ip and (|92p for the reduced nonequilibrium part of pair dis- 
tribution function A/ - j and the reduced nonequilibrium part of ionic potential A^j, respectively, 
defined by 

are graphically depicted in the case of £ = 1 in Figs. 2-3 to give pictorial representations for the 
nonequilibrium parts of the ionic liquid structure and the mean ionic potential in the Brownian 
motion model. They vanish as x and r increase to infinity from a finite value at the origin. The 
choice of the value of the reduced field strength £ is arbitrary; it could be as large as desired. 

Fig. 2 displays an important feature most distinguishable from the equilibrium pair distribution 
function for ion pair (J, i) that should be spherically symmetric and peaked at the coordinate origin 
(x,r) = (0,0). Instead, the nonequilibrium part of the pair distribution function A/^ (x, r, £) at 
£ > has a peak displaced from the coordinate origin. This means that the spherical symmetry 
originally present at equilibrium (i.e., at £ = 0) not only has been destroyed with its peak position 
displaced to a point (x, r) ^ (0, 0) from the coordinate origin, but also the ion atmosphere is no 
longer spherically symmetric if £ > 0. This means that the center of ion atmosphere has also 
been displaced by the the external field along the x axis. As a matter of fact, the present exact 
solutions of the governing equations provide the details of the state of distortion of the spherical 
ionic atmosphere and its migration as £ increases from £ = 0. We will see in the next section how 
this mode of distortion in the ion atmosphere is further modified in a manner of feedback process 
by the hydrodynamic motion of medium induced by the motions of ions under the influence of the 
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external field. Fig. 3 for illustrates the molecular cause for the distortion of the spherical 
ionic atmosphere through the nonequilibrium change in the ionic potentials. 

If the series representations for the Bessel functions Iq (uJir) and Iq (uJ^r) are used, the integrals 
can be evaluated in terms of elementary functions of £, x, and r, but since these series converges 
slowly, such series representations would have a limited practical value for precise evaluation of 
integrals. Nevertheless, such representation might be of some use for some theoretical study. Eqs. 
(|9ip and (|92p contain the information on the nonequilibrium ionic liquid structure and the mean 
potentials for the ionic liquid subjected to the external electric field. They are new results for ionic 
solutions in an external field examined here. The distribution functions and ionic potentials could 
be made use of to develop a theory of transport processes in binary electrolyte solutions. In this 
sense, they would be potentially very useful, especially, for calculating transport coefficients of the 
ionic solution in the electric field. 



IV. HYDRODYNAMIC EQUATION AND FLOW PROFILES 

In the conventional ionic conductance experiments the flow velocity of the medium is usually 
not large. Therefore flow may be regarded as laminar. Under this condition the nonlinear inertial 
term can be neglected in the Navier-Stokes equation. Moreover, the liquid may be considered 
incompressible to a good approximation. Under these conditions the Navier-Stokes equation be- 



comes the Stokes equation 19h21[] for an incompressible fluid. We therefore use the Stokes equation 
to calculate the flow velocity of the medium around the moving ions pulled by the external field. 
It may be helpful to point out that the flow field generated would be schematically reminiscent of 
the flow field around a moving object submerged in a medium. 

We assume that there are no body-forces other than an applied electric field. However, because 
ions are strongly correlated by long-range Coulomb forces and also interacting with the applied 
external electric field, it is necessary to calculate the mean local electric field. For the purpose of 
calculating it we may use the solutions of the OF equations and the Poisson equations presented 
in the previous section. Therefore the mean local electric field is expected to depend on the spatial 
position and the external field strength £. 
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A. Local Electric Field 



Since the field is aligned along the x axis and the charge density is given by the Poisson equation, 
the local force due to the field X on charge density g is given by 

DX. 



47T 



-V>i (r) • 



(94) 



Since ipj (r) is given by 



00 , cos (ax) ,._ , r , ,. , 
da ^-L [(R + 1) K (A1/0) 







R? 



+ (1-R) K (X 2 p) - 2 (1 - # 2 ) K (A 3 p)] 



+ 



a 



da sin (ax) [K (X 1 p) + K (X 2 p) - 2ivT (A3/0)] , 



(95) 



we obtain the mean local body-force 



XejK? 



+ 



da cos (ax) 
(1 - fl) (A| 



(R + l)(Xl 



a 



2R 2 



or 



X ej K 3 £ 
2ir 2 



2i? 2 

(ia sin (ax) 



-K (X 2 p) 



'-K (\ lP ) 
[1 - R 2 ) (A 



a 



R 2 



-Kq (X 3 p) 



a (A? 



a 



ft 2 



#0 (Aip) 



+- 



O: 



(Al - a 2 ) 



i? 2 



Kq (X 2 p) 



, "(^3-" 2 ) 
' « 2 



K (X 3 p) 



(96) 



This expression shows that the external force ejX is dressed up by the long-range correlations 
between the ions interacting through Coulomb forces and the interaction of ions and ion atmosphere 
with the external field. The effects of long-range correlations are described by the governing 
equations, and their feedback effect manifests itself in the form of dressed external force. This 
aspect is an important characteristic of the present theory of ionic solutions not usually seen in 
theories of charge carrier mobilities and their transport processes in recent literatures 
It is convenient to write Eq. (|96p in a compact form to solve the Stokes equation: 

.3 3 



Xe^K 



F x = 2 ^ 2 ^ [Q cos (ax) K (X t p) + Si sin (ax) K (X t p)] 



(97) 



1=1 
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where integral operators C\ and Si are denned by 

( (1+R) $-^ for/ = l 



da < 



(1 - fl) 2 (J- a2 ) for/ = 2 
.(^K^) for/ = 2 



(98) 



5« 



.-/() 



f ga(Af- a ') 

— w — 



da < 



€«(Aj-a 2 ) 
J? 



for I = 1 
for I = 2 



(99) 



- 2g ^r 2 ) for/ = 3 



This mean local force F x is an input for the Stokes equation of the flow problem under consideration. 



B. Stokes Equation and its Equivalent Form 



n 



At an arbitrary Reynolds number the steady Navier-Stokes equation[19|] must be used: 



p V • Vv - r/ V 2 v-?? fe V (V • v) = -Vp + F, 



(100) 



where ~p is the fluid density, r/o is the shear viscosity of the electrolyte solution, is its bulk 
viscosity, p is the pressure, and F is the body (external) force density. For an incompressible fluid 
V • v = 0. For a fluid undergoing laminar flow of low Reynolds number (typically Re = O(10~ 6 ) 
at the field gradient of 1 kVolt/m in aqueous solution) the inertial term can be neglected. Thus 
the Navier-Stokes equation for velocity v becomes the Stokes equations for divergenceless flow 



-??oV 2 v = - 
V • v = 0. 



-Vp + F, 



(101) 
(102) 



Note that the presence of an external field makes the pressure nonuniform in space. As is well 
known, if curl of Eq. (|10ip is taken, the Vp term vanishes and Eq. (jlOip takes the form 



r/ V xVxVxv = VxF. 



(103) 



Since VxVxv = V(V-v) — V 2 v by vector algebra, the two equations (jlOip and (|102p may 
be combined into a single equation 



?7oV xVxv = -Vj) + F. 



For the present problem F =S X F X , where S x is the unit vector along the x axis. 



(104) 
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To solve Eq. f)103|) for v, we observe V • v = 0, which means that there exists an axial vector 
A such that v = V x A, where A must depend on position vector r and field vector X, both of 
which are ordinary vectors. Thus we may transform the solution v of Eq. (|104p into the form 

v = VxVxa + v°, (105) 

where v° is a constant satisfying the appropriate boundary conditions of the velocity. Since a — > 
and also v should vanish as |r| — > oo, it follows v° = 0. Thus we will set v° = henceforth. 

In the first step to formally solve Eq. (|104p . substitute Eq. f)105j) with v° = into Eq. (|104[) to 
obtain the equation 

% V x V x V x V x a = - Vp + F. (106) 
By the identities of vector algebra 

V x V x a= V(diva)-V 2 a, (107) 
V x V x V (V • a) = 0, (108) 



and 



V x V x (V 2 a) = V (V 2 div a) - V 2 (V 2 a) , (109) 



it follows that 



VxVxVxVxa = VxVxV(V-a)-VxVx (V 2 a) 

= -V x Vx (V 2 a) . (110) 

Upon using Eq. (|110p in Eq. f)106|) and substituting the result into Eq. (|104p . we obtain a 
fourth-order differential equation of vector a: 

7/ V 2 V 2 a-F = V (7/ V 2 diva-p) . (Ill) 

This equation is equivalent to Eq. (jlOip or the Stokes equations. Because the left and right hand 
sides of Eq. (jllip are of two different kinds of vectors the equation may be separated into two 
equations: 

%V 2 V 2 a = F, (112) 
P = Po + % V 2 div a. 
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The solution of the Stokes equations is now reduced to that of Eq. (|112p , a fourth-order differential 
equation with F given by the solutions of the OF equations and the Poisson equations — namely, 
the governing equations. In summary, we have for the velocity and pressure the expressions 



v = VxVxa + v° = VxVxa, 
P = Po + %V 2 (V • a) . 



(113) 
(114) 



Vector a is determined by solving Eq. (|112p in terms of the local force density given by Eq. (|96l) 
or (|97p . a compact abbreviation of the former. In Eq. (|6ip po is a homogeneous pressure uniform 
in space, that is, the equilibrium pressure consisting of the osmotic pressure of the solution. This 
equilibrium pressure must be either supplied phenomenologically by using thermodynamics or from 
the statistical mechanics of the electrolyte solution. Therefore, given the solution for vector a, both 
velocity and pressure can be determined from the Stokes equation. 

To solve Eq. fll 12H for a, substitute Eq. (|97p into the former, which then reads 



V 2 (V 2 a) 



Xe k 3 

-S x [Q cos (ax) K (Xip) + Si sin (ax) K (Xip)] , 



2vr 2 r/ 



(115) 



where the repeated index I means a sum over / = 1, 2, 3. Since Eq. (|115p suggests that V 2 a must 
be a linear combination of the Bessel functions in the right hand side of the equation, recalling 
Eqs. flU and (|50D we find 

Xe^K 



V 2 a 



2ir 2 r] 



8 X 



cos (ax) K (Xip) sin (ax) K (Xip) 
T9 o r &l 



A 2 - a 2 



A 2 -a 2 



+ S X A*, 



(116) 



where A* is the homogeneous solution obeying the equation 



V 2 (V 2 A*) = 



(117) 



with A* = S X A*. The solution A* must satisfy the boundary conditions at infinite p. Thus we 
choose 



XejK 
2vr 2 r/ 



cos (ax) Ko(ap) sin (ax) Ko(ap) 
w T9 o 1~ T2 



A 2 - a 2 



Xf - a 2 



(118) 



since this satisfies Eq. (|117p . Therefore we obtain the equation 

Xe-:K 



V z &= 



2vr 2 7?o 



S T x 



cos (ax) [Kq(Xip) - K (ap)\ sin (ax) [K (Xip) - K (ap)] \ 

Cl Af3^2 + * Y2—2 



(119) 
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Since the solution of this inhomogeneous second-order differential equation must be a linear com- 
bination of the Bessel functions making up the inhomogeneous term on the right, it is sought in 
the form 

Xe* 



cos (ax) sin (ax) 



{&i [K (X lP ) - K (ap)] + b 2 [K (/3 lP ) - K (ap)]} 



(120) 



where 61, b 2 , and /3/ are constants determined as follows: On inserting this expansion into Eq. 
(TTT9D we hnd 



[61 (Xf - a 2 ) - 1] K (X lP ) + b 2 ((3? - a 2 ) K (f3 lP ) + K (a P ) = 0. 



(121) 



The expansion coefficients b\ and b 2 and the parameter are determined below. Since Bessel 
functions K (Xi P ), K ((3i P ), and K (a P ) not only do not vanish everywhere in p, but also their 
arguments are arbitrary, we may choose 61 and b 2 such that 

1 



A? -a" 



and 



lim b 2 (ft - a 2 ) K ([3 lP ) = -K (a P ). 



(122) 



(123) 



Then Eq. ([HI]) is satisfied and hence Eq. (fT20"j) is a solution of Eq. ([1191) . CCE1 therefore 
implies 

1 



(A 2 



a- 



'■) 



(124) 



Finally, we obtain for the solution of Eq. (|119|) 



r,^—S x [Ci cos (ax) + Si sin (ax)] x 

j [K (Xi P ) - Ko^of))] a P Ki(a P ) 



(A? 



2 (A? 



(125) 



For the solution (|125p for a, we have used Eq. (|107p and the recurrence relations of the Bessel 



functions 
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33] 



^K (z) = -K x (z), 
j-K x (z) = -K (z) - -K X (z) , 
1 d z±-l)zK 1 (z) = -2K (z), 



(126) 



z dz dz 

as well as diva = da x /dx owing to the fact that F =S X F X and hence a p = ag = identically. 
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C. Fourier Transform Solution for the Axial Velocity 



It is now possible to obtain the Fourier transform solution for the axial component of the 
velocity. Since 

(V x V x a) x = V x (diva) - V 2 a x , 

by using the formulas for V x (div a) and V 2 a x it follows from Eq. (|113p the Fourier transform 
solution for the axial velocity component for all values of x and p: 

3 r \2 



v, (x, P - = y Ci cos (ax) ± WW ~ KM! + 

^ 2^ V0 k^ 1 1 J \ (A 2 -a 2 ) 2 2(A 2 -«2)j 

jfgj A . f A 2 [if (Azp) -if (ap)] apK x (ap) \ 

> bi sin (ax J < — k 7—= > . (127) 

2tt 2 %«^ V ; \ (A 2 - a 2 ) 2 (A 2 - a 2 ) J k ; 



Here we now have restored the summation sign over index I. For a more explicit expression the 
sum over I may be expanded. This formula does not exactly agree with Wilson's expression 25|] for 
the axial velocity because of some missing terms and typographical errors in his formula. 

The Fourier transform integrals in Eq. (|127p may be expressed by using the reduced variables 
defined in Eqs. ([871) and ((551) to cast them into as simple forms as possible. We will also define the 
reduced velocity 

v = (2V2w 2 rio/zeXKj v. (128) 
Then the axial velocity v x (x, p; £) is given by 

(x, p; = " 7=TV v g (x, r, g) , (129) 

2V27T^7/o 

where the reduced axial velocity is now given by components made up of cosine and sine Fourier 
transforms: 

v x . (x, r, = \k* + 1=K? - Kt + \vKt - -Ltf|. (130) 
Various components in Eq. (j!30p are defined by the Fourier transforms 

K? (x,r,0 = |°° d * (1 C ! S 2^2 ) H K ^ir) + u> 2 2 K (u; 2 r) - ^Mk^t)] , (131) 



d , , , / , isin (art) 
tff (x, r, = / rft- x 



1 - 2£ 2 t 2 ) 

wfKo(wir) ^2-^0(^2^ 



+ 7 , v 7 v - 2w^o(w3r 



1 + Vi-2e 2 t 2 1- vi-2c 2 * 2 



(132) 
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K% = / dt cos (xt)u)iK (rt), (133) 
Jo 

poo 

Kl= dt cos {xt)tKi{rt), (134) 
Jo 

poo 

K%= I dtt sin (xt)K (rt). (135) 
J o 

Integrals -Kg, and i£| can be evaluated analytically in closed algebraic forms. On the other 
hand, the Brownian motion part of the integrals (x,r, £) and Kf {x,r, £) can be evaluated 
by methods of contour integration in the region satisfying Condition (190)) required by the Jordan 
lemma {30! for the contour integrals. Outside the region, they are computed by using straightforward 
numerical integration methods employing a method of principal values. 

1. Evaluation of Integrals K\, Kg, and Kl 

All the integrals appearing in the expression for v x (x, r, £) do not appear simple at first glance. 
Presumably, for this reason Wilson evaluated the integrals for the case of x = r = only. However, 
the integrals K%, K£, and K\ are indeed amenable to analytic evaluations in closed form. We 
explicitly illustrate the method by using K% as an example. Other integrals can be evaluated 
similarly. 

a. K\ On substitution of the integral representation [32|] of the Bessel function K u {rt) of 
integer order 

poo / . 

K u {z) = J dse- zcoshs cosh(vs) (\axgz\ < ~) , (136) 
the integral K\ can be written as 

-1 poo poo _ -, 

K%=-\ dt ds(l + t 2 ) e-Krccehe-ta) + e -i(rcoshs+ix) _ 

2 Jo Jo L ^ 

It is legitimate to interchange the order of integrals. Then the integration over t is trivial; changing 

variable to z = sinhs, we obtain elementary integrals with respect to z, which can be easily 

integrated: 

TT vr (2x 2 -r 2 ) 
K C A = ^ i -L. 137 

2(x 2 + r 2 ) 1/2 2(r 2 +x 2 ) 5/2 
It reminds us of Coulombic and dipole contributions, which are purely mechanical. 

b. K§ Upon using the integral representation of K\ (rt) and the same procedure as for integral 
-fT|, we obtain K£, 

ttt 

KZ = -—— j. (138) 
2(x 2 + r 2 )2 
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This integral appears in Wilson's formulation as the divergence-causing term. We will return to it 



again when we compare the present result with Wilson's 



251 ] in more detail. 



c. K\ This integral also can be evaluated in the same manner as for K%. We obtain 

K = 1?—,. (139) 

2(x 2 +r 2 )2 

The three integrals K%, K§, and K\ make up purely mechanical contributions to the axial velocity 
v x . They may be interpreted as either Coulombic or dipolar contributions of the ion atmosphere, 
which acts as if it is a dipole toward the external field. The collection of the three integrals 
evaluated up to this point will be collectively referred to as a mechanical velocity (y x ) , which is 
a countercurrent induced by Coulomb and dipole interactions of ion atmosphere interacting with 
the applied external field: 

(v..) mc EE -Kl + T -Kl - -LtfJ 



7T X 



2 v / 2£ (x 2 + r 2 ^ 



7T irr 



tt (2x 2 - r 2 ) 



19 + 3+— TF>' ( 140 ) 

2(2-2+^)1/2 4(x2+r2) | 2 {x 2 +r 2 f /2 

This contribution of (v a ) to v x represents the fully deterministic part of the hydrodynamic velocity 
that is not associated with the Brownian motion of particles giving rise to the dissipative part of 
the local body-force. In fact, one of these terms [i.e., the first term in the second equality of Eq. 
(|140p ], when inserted into the velocity formula (|129p . becomes field-independent and, consequently, 
does not contribute to the mobility or electrophoretic coefficient. Moreover, (v x ) me is negatively 
divergent at the coordinate origin, and its manner of divergence is clearly direction-dependent, 
that is, depending on whether the zero of x or r is approached first. Note that when converted to 
the axial velocity in real units, the last three terms in (v^) in Eq. (|140p are proportional to the 
reduced field strength £. 



2. Evaluation of {x, r, £) and Kf (x, r, £) Arising from Brownian Motions 

The remaining integrals (|13ip and (|132p can be calculated by means of contour integration 
methods described in Appendix A. We collect them in the form 

±Kf + -^Kf = ~ [C a (x, r, - 6t (x, r, 0] - \ [<t 2 (x, r, £) - 6 2 (x, r, 0] , (141) 
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where 



( i - y 2 + y/i + 2£V 



d (x, r; = / <*lP TTW -e _sew /o(E7ir), (142) 



o 



e>(.r.,-) - / dy * V ^I (oJ 3 r), (143) 

■ v^i+e 3 ) v^Cy (l - y 2 + v / lT2£V) 
e 1 (x,r;0= I dy ^ { e'^I^r), (144) 

(i + 2ev)(i + v / i+W 



o 

6 2 (x,r;£) = I* dy ^ V ^~fh -wi (u 3 r). (145) 



o 



1 + 2£ V 

The parameters (x, r) in integrals (|142p and (|144p are subject to the condition 



xy/l + C 2 > r (146) 

and integrals (|143p and (|145p to the condition 

x > r, (147) 

both of which arise from the condition to satisfy the Jordan lemmajsol on the contour integrals 
involving an infinite semicircle in the complex plane; see contours in Figs. 11-13 in Appendix A: 

x/r + Recoi (t)/lmt > (t = complex; I = 1, 2, 3) . (148) 

These conditions also have been mentioned in connection with the pair distribution functions and 
potentials in Sec. II. If these conditions are not met, the contour integration methods cannot be 
applied because the integrals along curve Coo of infinite radius diverge. In the region of (x, r) 
plane not satisfying these conditions (i.e., exterior to the region) the integrals must be evaluated 



numerically by using the method of principal values [341] . Note that as in the contour integration 
methods used for Eqs. (|142p - (jl45p the contributions from the singular points cancel in the end, 
leaving only the principal value parts. It is also important to note the sine transform terms make 
significant contributions comparable in magnitude to the cosine transform terms, as will be found 
later in the numerical analysis. On the other hand, if x were set equal to zero, the sine integral 
would identically vanish and thus have made no contribution to the velocity. Consequently, the 
final velocity values would be different depending on whether setting x and r equal to zero before 
or after integration. This subtle, but important point should be kept in mind when we handle this 
kind of integrals or the result obtained could be misleading. 
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In summary for the axial velocity, we obtain 
K 2 k B T x 



Vx {x, r, £) 



87T7/0 (x 2 + r 
zenX 



2^3/2 



1 



2x" 



(x 2 +r 2 ) 2 2(x 2 + r 2 ) 2 (x 2 + r 2 ) 
[Ci (x, r, + 2£ 2 (x, r, £) - 6i (x, r, - 26 2 (x, r, 0] 



(149) 



4-v/2vr?7o 

The axial velocity obtained here contains a term independent of the external field — i.e., the first 
term on the right, whereas the rest of terms are led by terms proportional to X (or £ in reduced 
units); they are in fact rather complicated functions of the field strength £. Physically, the velocity 
calculated from the Stokes equation represents the flow profile of the countercurrent induced by 
the moving center ion and its ion atmosphere in response to the external electric field. 



3. Electrophoretic Factor 

The mobility of ions in the x direction is associated with the field-dependent terms of the axial 
velocity, and the mobility or electrophoretic coefficient can be defined as the coefficient in the axial 
velocity vs. electric field according to the thermodynamic force-flux relations in thermodynamics 



35, 



36]. Therefore, according to the usual practice in the theory of ionic 



of irreversible processes 

conductance 0] within the framework of irreversible thermodynamics, we define I lie elect roplioret ic 
factor f (x, r; £) as follows: 

K 2 k B T x zeXn 

V* (x, r, £) = r - — = f (x, r; £)• (150) 

87r?7o (x 2 + r 2 )2 6V27T7/0 

In fact, the factor f(x, r; £) is generally dependent on position coordinates x and r as well as £. 
Note that the second term on the right of Eq. (|150f) is reminiscent of the velocity formula in Eq. 
(J3|) , which was obtained by a heuristic argument on the basis of the Stokes law in contrast to the 
hydrodynamic derivation of Eq. (|150p . Then upon comparison with the axial velocity formula 
p49D the electrophoretic factor f (x, r; £) is identified with the expression 

3 3r 2 3 (2x 2 - r 2 ) 

2(x 2 + r 2 ) 1/2 4(r 2 +x 2 ) 3/2 2 (r 2 + x 2 ) 5/2 

+ | [d (x, r, + 2£ 2 (x, r, £) - 6i (x, r, - 26 2 (x, r, £)] • (151) 

Since it generally depends on coordinates as does the axial velocity, the factor f(x,r;£), in fact, 
describes the electrophoretic profile in (x, r) plane as is evident from the figure shown below. 
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4- Numerical Evaluation of the Axial Velocity 

Since it is important to learn about the axial velocity profiles we have plotted them in the (x, r) 
plane in the case of £ = 1. In the region satisfying Conditions (|148p the formula given in Eq. (|149j) 
is used with the integrals (|142[) — fj!45|) . and in the exterior to the region defined by the conditions 
the velocity integrals for the Brownian motion contributions — i.e., Eqs. (|13ip and (|132p — are 
calculated by applying methods of principal integration because the integrals have singularities on 
the real axis. Thus computed axial velocity profiles are summarized in Figs. 4-6. 

In Fig. 4 the axial velocity is plotted in (x, r) plane in 3D with the vertical axis indicating the 
magnitude (color coded) of the axial velocity. It is seen negative in a semi-elliptic region enclosing 
the r axis beginning from r = (yellow-green color) as predicted by Formula (|149|) . it being negative 
principally because of the mechanical part of the axial velocity (v.j;) me , which becomes dominant 
over the Brownian motion contributions — the last group of terms in Eq. (I149p . According to Fig. 
4, the maximum of the axial velocity in the positive x direction (dark red region) is located in 
the neighborhood of the coordinate origin, but displaced from the origin (x,r) = (0,0). The axial 
velocity decreases gradually and eventually vanishes as x and r values increase to infinity. To gain 
a better idea of the electrophoretic factor f (x, r, £) it is plotted 3-dimensionally in (x, r) plane in 
Fig. 5 with the magnitude in the similar color coding used for Fig. 4. Its shape is rather similar to 
the velocity profile in Fig. 4, but its sign is opposite to that of v x owing to the way it is defined. To 
have a better idea of the behavior of the electrophoretic factor we have plotted the projection of the 
level curves of the f (x,r, £) surface onto the (x,r) plane in Fig. 6. It displays two sets of roughly 
elliptical contours, one with the major axis on the x axis and the other with the major axis on 
the r axis excluding the coordinate origin. The former set of contours corresponds to the negative 
portion of f (x, r, £) , whereas the latter to the positive portion of f (x, r, £) but transversal to the x 
axis. The outermost level curve denoted C p in fact represents the locus of zero of f (x, r, £), that is, 
f (x, r, £) = 0. These two sets of quasi-elliptical contours, and particularly, curve C p (i.e., the quasi- 
ellipse above the x axis) indicates how the spherical ion atmosphere at equilibrium with its center 
located at the coordinate origin when £ was equal to zero drifts away from the origin along the x 
axis and the spherical form is, at the same time, distorted to a non-spherical (quasi-elliptical) form 
with its center at x > as the external field strength increases — i.e., a nonequilibrium state. For 
example, in the present reduced variables employed, the equilibrium radius of the ion atmosphere 
(£ = 0) is (l/\/2) n~ l ~ 0.7/t _1 with the center at the coordinate origin, but if £ = 1, not only the 
center of the quasi-ellipse has migrated to x ~ 0.6k _1 and the curve C p is no longer spherical with 
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the major axis reduced to approximately 1.2k _1 = 2 (0.6k _1 ) instead of 2 (l/\/2j k^ 1 — IAk^ 1 in 
the case of £ = 0. This trend persists with increasing £. This behavior is numerically examined in 
Fig. 7, where the position of the center (x c , 0) of the quasi-ellipse is plotted against £. It gradually 
and significantly diminishes with increasing £ after having reached a maximum. Since this position 
(x c , 0) is at the center of displaced ion atmosphere that is simultaneously distorted by the external 
field it is natural to choose x in f (x, r, £) with the x coordinate of the center x c of the quasi-ellipse as 
the center of the ion atmosphere at £. Since the electrophoretic coefficient may be regarded as the 
force on the imaginary spherical ion atmosphere with its center at (x c ,0), then it is reasonable to 
choose r in f (x, r, £) with r = x c . With this choice of the x and r values in the electrophoretic factor 
f (x, r, £) we have verified that the electrophoretic coefficient thus calculated invariably produces 
the correctly behaved equivalent ionic conductance over a wide range of the external field strength, 
provided that the relaxation time coefficient [see Eqs. (|192p and (|193p below] is calculated with 
the same set of (x,r). Thus, in this manner we have been able to formulate a procedure based on 
computation result for selecting the position parameters (x, r) in the electrophoretic and relaxation 
time factors and therewith the ionic conductance unambiguously. We now state this procedure as 
follows: The values of the coordinates x and r in the electrophoretic factor f (x, r, £) are selected to 
be the x coordinate of the center of the quasi- elliptic level curve Cp and the corresponding value for 
r of the imaginary spherical ion atmosphere r c = x c centered at (x c ,0). The relaxation time factor 
is similarly calculated. In retrospect, this procedure — which may be called a rule — seems natural 
since the center of the ion atmosphere drifts along the x axis as £ increases and the electrophoretic 
coefficient must be reckoned with respect to the center of ion located at the (x c , 0) of the spherical 
ion atmosphere of radius x c , namely, (l/y/2) kT x in the actual units, which means r c = x c . 

With this identification of the coordinate parameters (x, r) in the electrophoretic and relax- 
ation time factors the electrophoretic and relaxation time coefficients are rendered unambiguous 
and unique. They are also divergence-free because the center of the displaced and distorted ion 
atmosphere does not occur at the coordinate origin for all values of £ and the OW theory becomes 
free from the divergence difficulty inherent to Wilson's procedure of selecting x = r = 0. 

5. Comparison with Wilson's Result for the Electrophoretic Coefficient 

Having defined the electrophoretic factor based on the full formula (|149|) for the axial velocity 
obtained from the Stokes equation, we investigate how Wilson's result for the electrophoretic 
coefficient can be recovered. He observed that since the ion of interest in conductance experiment 
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is the center ion of the ion atmosphere, which is located at the coordinate origin, the axial velocity 
must be considered at x = r = 0. He then noticed that the Fourier transform integrals comprising 
the axial velocity could be analytically evaluated at x = r = 0, because in the Bessel function 



Kq(z) represented in power series as 
1 



32, 



33] 



K (z) 



In | -z ] + 7 



+ 



2> ;2!r 



(152) 



, Q 2 2 {V.y 2 4 (2!)^ 2 6 (3!) 2 

where 7 is Euler's constant, if x = r = 0, only the leading term of Kq (z) contributes. Therefore, 
at x = r = Formula Q129)) for the axial velocity can be written as a sum of simple integrals 



u*(0,0;£) 



1 



1 



(1 - 2et 2 ) 

u\ In (^i) + ul In (^) - 2 (1 - i? 2 ) In (^) ] + efti. (153) 

The logarithmic integrals can be exactly evaluated by means of contour integrations by using 
contours similar to Figs. 10-12 in Appendix A. (Note, however, his contours used are not exactly 
the same as Figs. 10-12 we have employed in Appendix A except for the locations of simple poles 
and branch cuts.) With so evaluated integrals and the electrophoretic coefficient defined by the 
relation [3, 2jj] 



v x (0,0; 



zeXn 



no, 



6\/27rf?o' 

the electrophoretic coefficient /(£) could be shown given by the expression 

3 



(154) 



/(0 = i + 



4\/2£ 3 



2f sinh- 1 £ + V2£ - + e 



[1 + 2£ 2 ) tan" 1 (yfe) + (l + 2£ 2 ) tan" 1 [ J_ 



(155) 



provided that the last integral in Eq. Q153|) is ignored. For this formula for / (£) we have used the 
identities: 

X - tan" 1 fey/\ + = tan" 1 fe^/l + ( 2 ) \ sinh" 1 fey/ 1 + £ 2 ) = 2 sinh- 1 £. 

As a matter of fact, for the last integral in Eq. (|153p for v x (0, 0; £) Wilson[25] argued that the 
integral of \ contributes nothing to the electrophoretic coefficient because its contour integral 
vanishes. This argument is fallacious because although the contour integral in question certainly 
vanishes, it is composed of two integrals which are manifestly infinite, but opposite in sign: 

1 



-dz 



c 



00 L f 1 

—ax + / — dz 

00 2 Jc x 2 



0. 



(156) 
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As a matter of fact, according to the analysis leading to Eq. (|149p the last integral in Eq. (|153j) 
originates from the integral which we have already evaluated analytically for all values of x 
and r, and it is equal to zero at r = only if x 7^ 0, as is obvious from the following consideration: 

7rr as r —> 0, x 7^ 

Kt = w = { ^ . (157) 

2(x 2 + r 2 )5 for x -»• 0, r / 

However, if x and r simultaneously tend to at the same rate, K£ is manifestly divergent. There- 
fore rigorously speaking, Wilson's electrophoretic coefficient cannot be defined upon evaluation of 
v x (0, 0; £) with preset values of x = r = unless we simply abandon the divergent term. It now 
appears that his procedure of setting x = r = in the velocity integrals before evaluating the 
integrals is the cause for the divergence difficulty to obtain a finite electrophoretic coefficient, or 
the position x = r = should not have been taken in the electrophoretic coefficient defined through 
the thermodynamic force-flux relation for mobility or the Stokes law. This divergence difficulty 
and our desire to obtain physically sensible mobility coefficient was the principal motivation that 
we have evaluated and examined the velocity profiles in the (x, r) plane to understand how the 
velocity varies in space and to find out what would be the most probable or reasonable velocity 
that should be used to calculate ionic conductance if the Wilson-Onsager theory of conduction is 
adopted as the theory to rely on. We believe that OW theory is a correct approach to the ionic 
conduction problem, but the solutions must be evaluated more carefully for a wider range of (x, r), 
because the ion atmosphere, and therefore its center, migrates under the influence of an external 
field. 

We now would like to show in what manner the Wilson formula for / (£) would emerge from 
Eq. (|15ip . First of all, the potentially divergent mechanical contribution (v x ) me should be ignored 
to obtain a finite numerical value for f (x, r; £) at x = r = 0, although neglecting (v x ) me would 
result in a significant error to the axial velocity, and integrals £1 (x, r, £), • • • , 62 ( x , r , should be 
approximated as follows. First, let the Brownian motion contributions be denoted by 

VBrown (x, r, £) = £1 (x, r, £) + 2£ 2 (x, r, £) - 61 (x, r, £) - 26 2 (x, r, £) . 

Then if the Bessel functions Io(z) in the integrals for £1 (x,r, £), etc. are expanded in series 

00 (l z 2) k 
k=0 v ; 

they can be evaluated analytically in closed form (at least, for quite a few leading order terms) at 
x = 0. Especially, at x = the zeroth-order term, namely, the k = term in Eq. ()158ap . gives rise 
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to exactly the same / (£) as in Wilson's, Eq. ([1550 : 



VBrown (0, 0, £) = ~ fBrown(£)> ( 159 ) 



where f Brown (0 is then given by the expression 



3 /-v^OT) (i - y 2 + yi + gg) ,i 2cV(i-y 2 

f Brown (0 = / dy -, , oe2 2 H 3 / "2/ 



2 7 1 + 2^V y Q * 1 + 2CV 



3 f VW) V2fr l-y 2 + xATW /-i V2gy (1 - jg) 

2 -A) (l + 2^ 2 y 2 ) (1 + Vl + 2 eV) ^ 1 + 2£V 



which can be shown identical with /(£) in Eq. ([1550 . This process of arriving at /(£) from f(x, r; 
as a matter of fact, indicates that the electrophoretic coefficient /(£) obtained by Wilson[25] must be 
regarded as an approximation to the more precisely defined exact electrophoretic factor (or mobility 
coefficient) f(x, r; £) through the mobility 35j of ions on the basis of the irreversible thermodynamic 
force-flux relation between the external electric field and flow velocity. Recall that for /(£) it is 
necessary to leave out (va;) me from the axial velocity v x (x,r, £) in Eq. ([1491) . It is of course 
necessary also to leave out the field independent term — the first term on the right in Eq. ([1490 — 
for both /(£) and f(x, r; £) because the term has nothing to do with the mobility of ions in the 
external electric field. 

The axial velocity profiles presented in Eq. ([1490 arise from the presence of ion atmosphere 
and its interaction with the center ion itself and the external electric field. We must recognize that 
dynamics of ions in a solution and their interactions with the external field is not like that of an 
isolated single ion in the external field. Moreover, the center ion of the ion atmosphere does not 
directly contribute to the ionic conduction because of the countercurrent of the medium produced 
by the ion atmosphere, and the electrophoretic coefficient must be appropriately calculated taking 
this fact and the interaction of ion atmosphere with the external field into account. Therefore the 
position coordinate values should be taken with those of a point other than the coordinate origin, 
preferably, exterior to the curve C p , to calculate the electrophoretic coefficient because the center 
ion of the ion atmosphere moves with £ increasing; see Fig. 6 and the rule for choosing (x, r) in 
f(x,r;£) proposed. In this regard, recall that f(x, r;£) = on C p . 

D. Fourier Transform Solution for the Transversal Velocity 

By using the relation 

(V x V x a) p = V p (div a) - V 2 a p = V p (div a) (161) 
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in the case of a p = clq = we find the transversal velocity component in the form 



v ( x r , = -^=J? - -J? — 



J% + Jt — 2 r ^5 > 



where v p (x,r, £) is the reduced transversal velocity defined by Eq. (|128p . 

t 2 uJiKi{ujir) 
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00 cos (tx) 
o * (1-2^) 
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1 + v/l - 2£ 2 t 2 
t 2 u}\Ki{oj2r 
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dt , Sm : f L [twiA'i(wir) + tw 2 Ai(w 2 r) - 4^ 2 t 3 w 3 Ai(w 3 r)] , 



(i - 2^ 2 t 2 ) 



(162) 



(163) 
(164) 



dt cos (tx) tK\{tr), 
<it sin (te)t 2 Ki(rt), 



sin (tx) ti^o (rt) ■ 



(165) 
(166) 
(167) 



The integrals J|, J|, and J| are analytically evaluated by using the integral representations of the 



'4' U A 

Bessel functions in the same manner as for K%, and K^: 



■J i 



■I- 



nr 



2 (x 2 + r 2 )2 ' 
37rx (x 4 — x 2 r 2 — r 4 ) 
2r 3 (x 2 + r 2 ) 5 

7TX 



(168) 

(169) 
(170) 



2(x 2 +r 2 ) 3/2 ' 

The integrals and Jjf can be evaluated by using the contour integration methods similarly 
for the integrals and Kf . Jordan's lemma gives rise to the same conditions as Inequalities 
(|148p . In the region outside the validity of Ineq. (|148p the method of principal value integration is 
numerically employed. 

In summary, we obtain the transversal velocity component v p (x, r, £) in the form 

K 2 ksT r 



v p (x,r, £) 



8tTT]q ^ x 2 _|_ r 2 

KzeX 



+ 



4V27T?7 



xr 



+ 



3x (x 4 — x 2 r 2 — r 4 ) 



2(x 2 +r 2 ) 3/2 r 3( x 2 +r 2) 
[C 3 r, + ©3 (x, r, + 4£ 4 (x, r, - 46 4 (x, r, 0] , 



(171) 
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where 



U(x,r,t)= dy-* V 2 e-^hpir), (172) 

J o I 1 + z ? 2/ J 

€4(x,r,0= f dy ^Y^f -e-^h^r), (173) 

7o 1 + 2 ? y 



©3 (x,r,0 = / dy V V -e^I, for). (174) 

■/o ^2 (1 + 2£V) ^ + ^l + 2 eV) 

& 4 (x,r,0= dy V \ ^h^r). (175) 

These integrals can be analytically evaluated term by term by using the series expansion of the 
Bessel function I\(z) or numerically by a fairly straightforward procedure. The profiles of v p (x, r, £) 
look quite similar to those of the axial velocity v x . (x,r,£) in Figs. 4-6. 



E. Fourier Transform Solution for Pressure 

Since for the present system the (nonequilibrium) pressure is given by 

p = P0 + r]0 V 2 (^y (176) 

it is easy to calculate it from Eq. f|125[) : 

zeX f°° , a sin (ax) r/1 T ^ . 

+ (l-R) K (X 2P ) -2(l-R 2 ) K (\ 3 p)) 

_ zeX r 

zeX r°° a 2 cos (ax) 



4tt 2 



daa sin (ax) Ko(ap) 
£ da a2c °^ ax) [K (X lP ) + K (X 2P ) - 2K (X 3 p)} . (177) 



The formula presented above represents a nonequilibrium part of pressure Ap = p — po consistent 
with the velocity components obtained as the solution of the Stokes equation for a fluid in an 
external electric field. It also can be decomposed into the mechanical and Brownian motion parts 
as for the axial and transversal velocity components. They can be evaluated by the same methods 
as for the axial velocity, for example. With the reduced nonequilibrium pressure Ap defined by the 
formula 

. _ K ( zeXn 2 \ _l 

Ap = Api-^\ (178) 
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we obtain the nonequilibrium pressure profile in the form 



Ap 



TTX 



2 (x 2 + r 2 ) 3/2 



r VWW) <r xy y ( 1 + v^y + y/i + 2£V 

/ Wo, V 



') 



I (uir) 




(179) 



This shows that Ap is also singular at the origin of the coordinates. It is significant to observe 
that the nonequilibrium pressure Ap is generally negative, that is, there is a tension that becomes 
negative infinite at the origin. This implies that the nonequilibrium pressure is compressional in 
the neighborhood of the origin. Moreover, it is proportional to the field strength X. It seems to be 
a remarkable result, probably deserving a deeper consideration, because the degree of compression 
can be manipulated by the applied external electric field strength. We will report on a further 
study of this nonequilibrium pressure separately elsewhere [37)] . 

F. Ionic Field and Relaxation Time Effect 

Just as the velocity is induced by the mean local body-force which in turn is produced by 
interaction of the ion atmosphere [17] with the external field, the local ionic force field is modified 
by a feedback process of correlations arising from the Coulomb potentials and their interaction 
with the external field. Thus we may express the total electric field Xt acting on the ion in the x 
direction as 



where the local contribution AX is the ionic field produced by the interaction of the ion atmosphere 
with the external force field. If the potential of ion j in the electrolyte solution is denoted by tpj (r) 
the force arising from the potential ipj (r) is given by 



X t = X + AX, 



(180) 




(181) 



Upon using the solution for ipj (r), Eq. (|92p . we obtain the mean local ionic force 





(182) 
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in the notation already defined. It should be noted that the formula for ejAX (r) in Eq. (|182p is 
an exact result, although formal. The external field dependence £ enters the theory in a nonlinear 
manner through the arguments of the Bessel functions. We have shown that the Fourier transforms 
such as those in Eq. (|182p can be reduced to finite quadratures consisting of regular Bessel functions 
of zeroth order of second kind Iq(z) weighted by some algebraic functions. 

The integrals in the expression for AX (r) were also evaluated by Wilson in his dissertation 
for the case of x = r = in the same manner as for the electrophoretic coefficient /(£), Eq. (j!55|) . 
With the so-obtained result the relaxation time coefficient g (0 was defined by the relation 

ejfj/n 



zeAX (0, 0; f ) 



2D 



-9(0 



2D 



9(0, 



(183) 



for which he obtained g(0 in a simple analytic form 



9(0 = ^3 



£^0+0 - tan" 1 ( y=f 1 " ^ + tan' 1 (^) 



(184) 



It is a nonlinear but well-behaved function of the reduced field strength £; its limiting values are 
g(0) = | (2 — y/2) as £ — > and 5(00) = as £ — > 00. We will show presently under what condition 
this result, Eq. (1184p . is recovered from the exact formula for ejAX (r), Eq. (|182j) . 

To obtain a complete formula for the local ionic field from Eq. fll82|) the Fourier transform 
integrals therein must be calculated without setting x = r = before evaluating them. For this 
purpose we use the same contour integration methods described in Appendix A as for the velocity 
formulas. We thereby obtain ejAX (r) in the form 



ej AX (r) 



ejK 2 ^ 
2V2D 



V 2 (!+? 2 ) e-*V/ far) 
V 1 + 2£V 



2 / dy 
10 



e xy y 2 Io(uj; 



u 3 r) 



1 + 2£ V 



2 2 
AD 



yff&F) e-^y (l + y/1 + 2ey 2 ) k far) 

dy 1 i o£2 2 

1 + 2£ V 



+4£ 2 



1 dy y3e ' XVI ° (ZJ3r) 



(185) 



/o 1 + 2£V 

The variables x and r as well as £ and y appearing in the integrals in Eq. (|185j) are dimensionless 
reduced variables defined earlier. The integrals in Eq. (|185D . of course, are subject to conditions 



deduced from conditions (|148p related to the Jordan lemma [301] for the contour integrals. Exterior 
to the region of (x,r) satisfying Conditions (|148p the method of principal values for singular 
integrals is used to numerically compute the Fourier transform integrals. 

The first group of terms in Eq. (1185P descends from the cosine transform terms in Eq. (I185P 
whereas the second group originates from the sine transform terms. If x and r are set equal to 
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zero in Eq. (|185p the first group exactly gives rise to Wilson's result, Eq. (|184p . but the integrals 
in the second group do not vanish even if x = r = is taken after their evaluation, but give rise to 
a field-independent term —ejK 2 /2D in the limits of x — > and r — > 0, in addition to ^-dependent 
terms as shown below. However, if we took x = in the sine transform integrals in Eq. (|182p there 
would have been no contribution from it at all. This example, once again, manifestly demonstrates 
a need for caution to take in evaluating the Fourier transforms, especially, with regard to setting 
x = r = 0: We reiterate that the values obtained of the integrals are different, depending on 
whether particular parameter values, especially, x = r = 0, are taken before or after evaluation of 
the integrals, or even depending on the order of taking the limits x — > and r — > 0; we have seen 
a similar situation in the previous section for velocity profiles. The fact that x and r must be set 
equal to zero to obtain Wilson's formula for g(£) also suggests that his formula is an approximation 
to the full relaxation time factor defined below, as is /(£) an approximation to the full f (x, r;£) 
given in Eq. (|15ip . 

To define the appropriate relaxation time factor we split ejAX (r) in Eq. (|185p into two parts 
ej AX (r) = - ej ^-9c (x, r; - ^9s(x, r; (186) 



with the definitions 



1 /V2(i+a e -V/o (wir) /- f 1 e- x y y 2 I for ) 



1 ry/m?) e-wy f 1 + VT+2fV) (wir) 
9s(x,r;£) = - dy 



1 + 2£ V 
,2 f 1 y 3 e- x yi {u; 3 r) 

The reason for the splitting made above is that the term g s (x,r;^) tends to a value independent 
of £ as £ -)■ 0: 

5s (x,r;0) = ^ dye- xy yl (rs/2 - y 2 ] ^ 0, (189) 
which contributes a field-independent term to ejAX (r): 

J^dye- x yyI (rV2-y 2 ), (190) 



2D 

and this contribution would have nothing to do with the relaxation of ion atmosphere. Therefore 
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it is useful to define the field-dependent part of g s (x,r;£) by the expression 



Ag s (x,r;£) = - 



1 rVW^e) e~ x y y (1 + y/1 + 2£ V) J (SJi7 

2 7 dy " TTW 



dye iy yIo [ry/2- y 2 



y 3 e xy I (uJ 3 r) 

+ 2 £ / dy — i iota 2 
'o 1 + 2 S V 



With this we are now able to cast e,-AX (r) into a more appropriate form 



ej AX (r) 



2D 



g s (x,r;0) - J 2p-8(x,r;£), 



(191) 



(192) 



where the relaxation time factor g (x, r; £) in the cylindrical coordinate representation is defined 
by the formula 



(x, r; = y c (x, r; £) + A^ s (x, r; £). 



(193) 



It is easily verifiable that Ag s (x,r;£) is indeed a constant in the limit of £ = and hence 
£Ag s (x, r; £) vanishes as £ — >• 0. With this definition of relaxation time factor q(x, r;£), we are 
now ready to examine its relation to Wilson's <?(£) formula. 

If both x and r are set equal to zero in integrals in Eqs. (|187h . (I190p . and (|19ip . it follows 

dy- — ^ 

'o 



9c (0,0; 



1 



/V 2 ^ 2 ) y 2 ,1 

7o dy TTW?-^ 2 L dy T 



^2j ""1 + 2^V ./o a l + 2£V 

^3 Vf 2 + 1 - arctan 2f \/C 2 + 1 - 2 ^£ + 2 arctan ^) . 



(194) 



A<? s (0,0;£) 



1 / rV 2 ^) y 1 + vT+Wy 2 ) \ r 1 



4^3 ( t ln ( 2 ^ + 1) + ( 2 ^ + 1) " 1] " ^) + 4^3 t 2 ^ " ln + !)] 
0, 



5s (0,0;0) = 1. 



(195) 
(196) 



Here g c (0, 0; £) is identical with Formula g(£) in Eq. (|184p for the relaxation time coefficient 



in Wilson's method 



251 ] . whereas Ag s (0,0;£) together with g s (0, 0; 0) represents extra terms not 



present in his result. We reiterate that in Wilson's method g s (0,0;£) does not appear because the 
sin (ax) term in the sine transform integral vanishes if x is set equal to zero before evaluating the 
sine transform integral. This shows under what condition Wilson's g(£) is recoverable from the 
present result for g (x, r; £). 
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In Fig. 8, g (x,r;£) computed, for example, at x = 0.5 and r = 0.5 is plotted as a function of 
the field strength £ and compared with Wilson's relaxation time coefficient g (£) in Eq. (|184[) . the 
dotted curve. To better comprehend the profile of the ionic field graphically, we plot a 3D example 
of g (x, r; £) in the case of £ = 1 in Fig. 9. A combination of formulas ([193p - (|196p and the method 
of principal values for integration is used to compute the relaxation time factor g (x, r; £) presented 
in Fig. 9. As does the electrophoretic factor, it also exhibits a singular behavior near the origin, 
although the details are different from the behavior of f (x, r; £) in Fig. 5. It also vanishes as x and 
r increase to infinity. 

We give a short summary of this long section: The results of evaluation of v x (x, r, £), v p (x, r, £), 
and Ap (x, r, £) for all values of x and r and the related electrophoretic and relaxation time factors 
constitute some of the important contributions of this work to the hydrodynamics of strong binary 
electrolyte solutions in the external electric field. On extensively studying the velocity profiles we 
have been able to formulate a rule for selecting the position variable (x, r) in the electrophoretic 
and relaxation time coefficients, which are finite everywhere. By using this rule and the profiles 
of velocities and nonequilibrium pressure as well as the distribution functions and mean potentials 
calculated, we will also be able to predict or deduce, in a well-defined manner, hydrodynamic 
consequences to transport properties, such as conductivity, and related nonequilibrium properties 
of ionic motions in the medium in an external electric field of arbitrary strength. 



V. DISCUSSION AND CONCLUDING REMARKS 



In this paper, we have shown that since ions interact with each other through long-range Coulom- 
bic interactions, ion atmosphere with ions, and both of them with the external electric field, the 
correlations of particles in the ionic liquids are quite complex and the whole body of an ionic so- 
lution collectively and cooperatively moves subjected to the external electric field. Consequently, 
even at a dilute ionic concentration the macroscopic behavior of electrolyte solutions under an 
external electric field is not simple, but rather complex and, therefore, exhibits an interesting feed- 
back system. In this regard, the subject matter is interesting from not only the theoretical, but 
also practical standpoint to gain insights into the behavior of complex liquids. For this reason, we 
believe the ideas of the OW theory [16;, 



22, 



251 ] as a theory designed to treat ionic fluid systems in 
the external field are worth studying in depth for the insights they provide for dynamical theories 
of ionic matter in general. However, examining in detail the solutions of the Stokes equation for 
flow velocity obtained from the solutions of the governing equations in the OW theory, we find that 



43 



the velocity formula not only had a divergence difficulty that we have unexpectedly encountered 
while studying it, but also was incompletely treated mathematically in Wilson's work[25] because 
only the behavior of the center ion at the coordinate origin was examined despite the fact that the 
ion atmosphere moves in the external electric field and develops a non-simple spatial structure. 
Therefore, we felt that there still remained the task of fully implementing the theory in a math- 
ematically satisfactory manner to make it serve as a complete theory of ion conductivity in the 
nonlinear regime of external field dependence. 

To achieve the goal in mind, we have numerically studied the velocity profiles in the configuration 
space over a range of external field strength and, in particular, the movement and distortion of 
the ion atmosphere, as the external electric field strength is continuously varied over a wide range. 
Thereby we have numerically quantified the trajectory of the center ion of the ion atmosphere with 
respect to £, but also studied the manner of its distortion from a spherical form to a quasi-elliptical 
form, as the field strength £ is varied. The general picture we obtain of the electrophoretic factor 
f (x, r;£) is as follows: within curve Cp it is negative whereas outside Cp it is positive. Moreover, 
Cp is non-spherical. This implies that the ion atmosphere not only polarizes into a negative and 
a positive domain (typical of a dipolar distribution), but also the boundary curve (i.e., Cp) gets 
distorted to a quasi-ellipse from a spherically symmetric form, as £ increases from zero. On the 
basis of the body of numerical studies of the axial velocity profiles we have been able to formulate 
a procedure by which it is sufficient to calculate the center position of the moving ion atmosphere 
at every value of £ and therewith calculate the electrophoretic and relaxation time coefficients as 
functions of £. 

The identified procedure is that: the electrophoretic coefficient at a value of £ is given by the 
electrophoretic factor f (x,r;£) evaluated at the center of the displaced spherical ion atmosphere of 
radius x c , whose center is located at (x,r) = (x c ,0), the center of the displaced quasi- elliptic curve 
Cp that is the locus of f(x, r;£) = 0. Since the spherical ion atmosphere with its center at (x c ,0) 
has a radius x c , the value of r c must be equal to x c . Therefore the electrophoretic coefficient f (£) 
is given by f (£) = f (x c ,r c ;£) = f (x c ,x c ;£) according to this finding. Since the center position of 
quasi- elliptic curve Cp is unique for every £ ; the electrophoretic coefficient f (£) defined is unique. 
The relaxation time coefficient (£) is then calculated by q (£) = 0(x c ,r c ;£) = fj(x c ,x c ;£) to be 
consistent with the electrophoretic coefficient defined. 

This behavior (trajectory) of the center of ion atmosphere gives rise to non-divergent elec- 
trophoretic coefficients for all field strengths and hence the ionic conductance based on the Fokker- 
Planck equations employed is now rendered divergence-free. This is made possible by recognizing 



44 



that the electrophoretic coefficient must be calculated for the moving ionic atmosphere with the 
center of the displaced ion atmosphere at (x c , 0) when the external electric field is applied to the 
system. 

The set of values for x and r obtained to use for f(x, r;£) and g(x,r;£) is (x c ,r c ) = (x c ,x c ) 
which is in significant contrast to the values x = r = taken by Wilson to evaluate the integrals 
that gives rise to a divergence difficulty. In the companion paper [31[, we apply this identification 
of x and r to compute the electrophoretic and relaxation time coefficients, and calculate therewith 
the equivalent ionic conductance and, in particular, the Wien effect of a binary electrolyte solution 
in comparison with experimental data available. 

The velocity profiles graphically presented also suggest a skin effect by which the mobility of 
ions in solution is predominantly contributed by ions outside the curve Cp. We have not suspected 
the existence of this effect before: that the conduction currents are mostly carried by ions in the 
periphery — i.e., in the shells of radius of O — of ion atmosphere, but not by the center ions, 
as is obvious from Fig. 4-Fig. 6. 

The another important mathematical question we are answering in the present work is that 
variable parameters, such as x, r, and £, in the Fourier transform solutions of the OF equations, 
Poisson equations, and Stokes equation should not be set equal to zero before fully evaluating 
them, since the results so obtained do not generally yield the same results as those obtained by 
setting them equal to zero after their complete evaluation. They would give identical results only 
if the results of the integrals are analytic everywhere in the space of x, r, and £, but the examples 
we have studied definitely show that the evaluated results are not necessarily analytic everywhere 
in the aforementioned space, and as a consequence the results of evaluations by the aforementioned 
two different modes can be significantly different; that is, the results are not uniformly convergent 
to the same conclusion. This should be regarded as a significant point of the present analysis to 
keep in mind in the study of this line of theories for ionic solutions. 

What we have shown in this work are the exact velocity and pressure profiles in space in a 
Brownian motion model, which we may apply to study other irreversible phenomena in the binary 
electrolyte solutions in the electric field than the Wien effect. Being full exact solutions without 
an approximation within the framework of the Brownian motion model, not only do they, at least 
in the low density regime, promise to provide a more complete picture of conduction phenomena, 
but also the insights gained therefrom should also help us develop theories of related trans por t 
phenomena in systems [S- 13] of current interest in science and engineering, such as plasmas^, 



semiconductors [4i. la. l39l] . small systems [40[|, etc. in electromagnetic fields. In any case, they repre- 



45 



sent new results in the subject field. In the sequels [14l] we will also study asymmetric electrolyte 
solutions, in which charges of the cation and anion are not symmetric, and ionic conductance 
under an external electric field. Although more complicated than the present symmetric binary 
electrolyte solutions, we find that a similar mathematical analysis is possible to obtain for them. 
The results of the mathematical solutions will be reported in the near future [l^. together with 
their numerical results [15] in comparison with experimental data. 
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Appendix A: Contour Integration Method 

Some of the integrals in the formulas for the distribution functions, potentials, velocities, and 
pressure can be evaluated by applying the method of contour integration which yields formulas 
more readily amenable to analysis and further approximations giving rise to simple results which 
will make them possible to use for the purpose of assessing the existing results on the subject 
matter. 



1. Axial Velocity 

We consider the axial velocity first for the reason that it contains more experimentally direct 
features than the nonequilibrium structure and potentials. The integrals appearing in the formal 
Fourier transform solutions in the present theory all involve the Bessel functions Kq(Xip) (I = 1, 2, 3) 
of argument Xip with A; being relatively complicated functions of the integration variable, the wave 
number; see Eq. (|87p and Eq. (|88p and also Eq. (j89|) for tJ;. In reduced variables we have defined 
for the analysis they have the mathematical properties listed below. 

(1) The zeros of the arguments of the Bessel function K u (uir) for r/0 are found to be : 

t = ±V2(1 + £ 2 ) for wi, (Al) 

t = for u 2 , (A2) 

t = ±i-^= for «3. (A3) 
v2 

The argument of K v {uj\r) therefore has branch points at tr = ±i-y/2 (1 + £ 2 )r, whereas the 
argument of K u {u)2 r ) has branch points at tr = x r and — oo x r and the argument of 
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K v {uj^r) branch points at tr = ±ir. 

Thus we may insert a branch cut on the imaginary axis of complex t plane between t = 
iy/2 (1 + £ 2 ) and t = —i\j2 (1 + £ 2 ) for the integral of K u (uixr), while a branch cut may be 
inserted along the negative real axis for the integral of K v (ui2r), and on the imaginary axis 
between t = i and t = —i for the integral of K v (uj^r), respectively. See Figs 11-13 below. 

(2) We recall that Bessel function K v (z) of complex variable z is regular in z plane cut 



along the negative real axis 



32, 



331 ] . That is, it is a multi- valued function in the cut 
plane. Therefore, in the present case, Ko(ojir) changes discontinuously as the branch cut 
— iry / 2 (1 + ^ 2 ), ir\j2 (1 + £ 2 ) is crossed (r is a fixed parameter), whereas K v (u<iv) changes 
discontinuously as the negative real axis is crossed, and K u (oj^r) changes discontinuously as 
the branch cut [—ir, ir] is crossed on the imaginary axis of t plane. Note that the Bessel 
functions K v (oj2r) and K v (tr) are defined in t plane cut along the negative real axis. 

(3) We also observe that all the integrands of the singular integrals, for example, in Kj? and Kf 
in Eqs. (|13ip and f|132|) are even with respect to t. 

(4) Moreover, for < argt < ir we find 

-(|-*) <argw,< (!>*>0; Z = 1,2,3). (A4) 

Therefore in the upper half of complex t plane 



lim K v (coir)= lim J— e -w,r -»• 0. (A5) 
|t|-»oc |t|-»oo V 2uir 

(5) Lastly, all the integrands in Eqs. (|13ip - (|132p have simple poles on the real axis at 

t = ± _L m 

There is also a branch cut between t = — and t = because of the (1 ± R) factor in 
the integrals, but this particular branch cut associated with \J\ — 2£ 2 i 2 does not play a role 
in the contour integrals considered in the present work, because the real axis is not crossed 
by the contours in performing integrations. Therefore we may ignore this particular branch 
cut. 



All these properties (l)-(5) to 



methods of contour integration 



aether suggest it is possible to evaluate the integrals by using 



301 ] along the closed contours of a semicircle as depicted in Figs. 
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8-10. However, in this approach the results obtained would not cover the entire region of the upper 
positive quadrant of plane (x, r). In the region outside the domain defined by the inequalities, Ineq. 
(|148p the Fourier transform integrals must be computed numerically because it is the region where 
Jordan's lemma 3(| is violated; that is, the contour integral along the circle does not vanish. 
In the exterior region their numerical values are small and hence of no importance. The practical 
advantage of this kind of contour integration method of evaluation is to isolate out the major part 
of contributions to the integrals and discuss the connection with the existing results where possible 
and with experimental data. It would be convenient to decompose Kf and Kf into component 
integrals as follows: 

Kf (x, r, = Kf + Kl - 4£ 2 if 3 c , (A7) 
Kf{x,r,i)=Ki + Ki-ZKl (A8) 

where Kf and Kf are in the order of their appearance in Eqs. (1131 1) and (|132p . 

Since methods of integration will be similar for the integrals involved in Kf and Kf (I = 1, 2, 3) 
we will illustrate them with the examples of integrals in Kf and Kf in the following. The results 
for the rest of integrals can be similarly obtained. 



a. Contour integrations of Kf and Kf 



As prototypes of contour integrals appearing in the axial velocity formula, integrals Kf and Kf 
are explicitly evaluated below; see Eq. ()130p -Eq. ()135p . Integrals Kf and Kf both have simple 
poles at t = ± (V2£)~ . There is a branch cut along the imaginary axis between t = -1^2(1 + £ 2 ) 
and +i\/2 (1 + £ 2 ) and also a branch cut on the real axis between t = — 1/\/2£ and t = +l/\/2£, 
but the latter branch cut plays no role in integration since the path of integration stays above the 
cut. For this reason the latter branch cut is not shown in Figs. 11-13. For evaluation of both Kf 
and Kf the contour in Fig. 12 is used. 

Consider the contour integral denoted by C\Kf along the contour C\ in complex plane z depicted 
in Fig. 9: 

dKf= I dze ix \ "} c2 2 K (Lj ir ), (A9) 

where 



1/2 



1 + z 2 + y/i - 2£ 2 z 2 . (AW) 



48 



Since there is no singularity enclosed by contour C\, this contour integral C\K\ is clearly equal to 
zero. Integral C\K^ can be decomposed into integrals along the paths C_, C+, C, Coo, and along 
the real axis t in Ci. We thus may write it as 



,00 1 + t 2 + v/l - 2£ 2 t 2 
C!jAT ~ ' J U 1^2^ L( -* tKc (l)r) 



l + z 2 + s/l - 2£ 2 ; 



C- 



1 + z 2 + v/i - 2e 2 ^ 2 

1 _ 2 g 2z2 L e lxz KviuJx (z) r) 



c 



1+Z 2 + y/1 - 2i 2 Z 2 

dz ~ 1 - 2i 2 z 2 Le%XZR ^ ( z ) r ) 



1 + Z 2 + y/l - 2i 2 Z 2 

+ / L ^ Z K^ (z) r) 

J Coo ^ 

= 0. (All) 
The first integral on the right, the integral along the real axis, can be shown to be equal to 2K^: 



OG 



l+t 2 + Vl-2£ 2 * 2 , 
x dt ~ 1^2^ ^Mui (t) r) 

/°° ui 2 
dtCOS {m) (1 - 2\ 2 t 2 ) K °^ ir) 



2K\. 



The remaining integrals along contours C_, C+, C, will be denoted by C-Kf, C+iff, Ciff, 



C Q K^, respectively. By the theorem of residues 30| the integral C_K^ gives tri times the residue of 
C-K{ at t = -1/V2& 



C-KZ = i^- (2f + 1) e-/^ Ko( - r) 



where 



Similarly, we obtain 



00 = — = — . 
V2i 



(A12) 



C+Kf = -m^ (2f + 1) ^/^KoOStr) 
Thus combining the results for C-K\ and C+if^, we obtain 



V2tt (1 + 2£ 2 ) 



.r 



C-Kt + C + i^ c = ^ 3 y sin ^— ) K (57r) . ( A 1 :! ) 
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To transform the contour integral CK\ around the branch cut along the imaginary axis we 
observe that if the phase of the argument of Kq (oo\r) on the right hand lip of the cut is chosen 
equal to zero, the phase of the argument on the left hand lip is %i, so that the argument has the 
form e m u)\r for the Bessel function on the left side of contour C. In this connection, it must be 
recalled that only the relative phase across the branch cut is of importance. When traced along C 
from the left to the right side of the cut, the Bessel function must be continued from the left side 
of the cut to the right side by the following continuation formula 



32 



33] 



K (e™z) = K (z) - TTiI (z) , 



(A14) 



where Jo (z) is the regular solution for the second kind of the Bessel function of order 0; Kq(z) is 
irregular in contrast to Jo (z). The irregular Bessel function K u {z) {y > 0) diverges logarithmically 
as z — > 0, whereas in series representation the Bessel function Iq (z) is regular and given by the 



formula 



32 



33] 



*>(*)=£ 



1 \ 2m 



m=0 



(m\) 2 



(A15) 



This function is finite at z = 0, but it behaves asymptotically as 

I Q (z) ~ (2irzy 1/2 e z [l + O^" 1 )] (\axgz\ < ^ttV (A16) 
Using formula (|A14|) and changing variable from iy to y, we obtain the integral along contour C: 



V5o+?) e~ x y (i-y 2 + ^/T+Wy 2 ~ 

( r: I dy V i - ^ 2y2 '-In r ) 



where 



l-y 2 + vr+2ev 



1/2 



(A17) 
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If the series form for Iq(z) in Eq. (IA15D is used, CKf can be computed in terms of quadratures 
of elementary functions — in fact, incomplete Laplace transforms. It should be recalled that this 
integral (|A17p is subject to the condition (|148p for the relation of x to r that is deducible from the 
Jordan lemma jsol]. To satisfy this lemma the integrands of integrals in K*? and Kf must satisfy 
the condition 



x Im t + r Re uj^ > 



1,2,3) 



(A19) 



for x, r > 0. Thus values of x and r in the (x, r) are limited to the region satisfying this condition 
plane, assuring the contour integrals along the curve Coo vanishes as x and r tend to infinity. In 
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the case of integral Kf the condition implies the inequality 



(A20) 



which in fact assures that the integral vanishes as x and r tend to infinity. Outside this region 
the contour integration method is not applicable. Therefore, integral (IA18j) does not hold and the 
Fourier transform integral Kf must be evaluated numerically. However, the numerical values of 
the integral in the exterior region gets diminishingly smaller as x and r increase to infinity. The 
contour integral along the outer semicircle does vanish in the region satisfying Jordan's lemma. 

Collecting the results for the contour integrals obtained above into Eq. (|Alip . we obtain the 
integral Kf in the form: 



Kf 
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8£ 3 

~ 2 7o 

The procedure of evaluating integrals Kf with the contour in Fig. 13 is entirely parallel to the 
one presented above for Kf. The result for the integral Kf is 

7T (1 + 2£ 2 ) 



Kf 



8£ 4 



cos (xt) Ko(ur) 



7T 
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This integral also is subject to condition (|A20|) and in the exterior region the Fourier transform 
integral must be evaluated numerically. Therefore, collecting results for Kf and Kf, we obtain 



-Kf + 4=^f 
2 1 V2 1 



7T(l + 2g 2 ) 

V / 20+a 





sm 



V2i 



cos (xt) Ko(ujr) 



K (ZJr) 



e-*y [l-y 2 + Vl + 2£V 
dy : _ - „ — x 



1 



1 + 2£ 2 y 2 
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(A23) 



i + yrr^v^ 

The first line involving sine and cosine functions in this result will be seen canceled by similar 
terms in the results for integrals K^ K^ K%, and 



b. K% and K s 2 



Evaluation of these integrals proceeds similarly to that of Kf and Kf with the contour given in 
Fig. 12 except that since the integrand does not have a branch cut on the imaginary axis, there 
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is no integral along the imaginary axis. There are only contributions from the residues at the 
singularities. They give rise to the following results: 



K c 2 
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Therefore we obtain 
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K (UJr). 



(A26) 



c. and J^f 



In the present cases, the integrands involve a branch cut along the imaginary axis from z = — i 
to +i. Therefore the appropriate contour to use is C3 depicted in Fig. 13. Evaluation of integrals 
-K3 and -K3 is entirely parallel to those of integrals K\ and Kf. The results of their evaluation are 
as follows: 



Ld 3 r) 
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2 1 y— TT^V - ' 



Therefore we find 
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(A29) 



2 7 - 1 + 2CV 

The integrals in Eqs. (|A27p and (|A28j) are subject to the condition deduced from the Jordan 
lemma, namely, x > r. 



d. Summary for the Reduced Axial Velocity 



Collecting the results presented earlier, we obtain the reduced axial velocity 



TT 



j [Cj (x, r, - 61 (x, r, 0] - I [£2 (x, r, £) - 6 2 (x, r, 0] , (A30) 
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where various symbols are denned by 





(A31) 



mc 
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As noted earlier, the terms made up of trigonometric functions in Eqs. (|A23j) . (jA26]), and ()A29p 
indeed cancel each other out. This velocity formula (IA30j) is the axial velocity profile of the 
countercurrent of the ion and its ion atmosphere in the coordinate system fixed at the center ion 
of the ion atmosphere, both of which are pulled by the external electric field. The first four terms 
making up (v x ) me (x,r,£) represent the "deterministic" part of the velocity v x , and the integrals 
d, £2, ©i, and ©2 involving the Bessel functions Io(uj\r) and Jo(^3 r ) stem from the Brownian 
motion part of the mean local force — namely, the dressed-up part of the local body force arising 
from the interaction of the center ion, its ion atmosphere, and their interaction with the external 
electric field, which distorts the ion atmosphere to an asymmetric form. This velocity formula 
obtained in Eq. (|A30|) is i n a convenient form to analyze Wilson's result, further examine the 
cause of divergence, and find a way to avoid the divergence difficulty in the evaluation of the 
electrophoretic coefficient. This aspect is discussed in the main text. 

2. Distribution Functions and Potentials 

The same methods of contour integration can be employed for the distribution functions /ji(r) 
representing the nonequilibrium ionic liquid structure (pair distribution function) and the mean 
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ionic potential ipj (r). They are summarized below: 
fij = fji (± r ) 



and 
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+ 2£V 

in the region of (x, r) in the upper plane where Jordan's lemma is satisfied. These results can be 
easily obtained by using the contour integration method described earlier in this Appendix. We 
notice that the nonequilibrium pair distribution functions and potentials do not contain mechanical 
contributions, but only the Brownian motion contributions. The reason is that fji{r) and ipj (r) 
are solutions of the OF equations and Poisson equations for the nonequilibrium part described by 
the Brownian motion model. 
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Figure Captions 



Fig. 1 The cylindrical coordinate system employed. The x axis is parallel to the external 
electric field. 

Fig. 2 Nonequilibrium part of the distribution function A/^ (+r) is plotted in (x, r) plane 
at £ = 1. Here Aj^ (+r) = [fa (+r) - n 2 ] / (V2nze 2 /irDk B T) . Af {j (+r) is computed with the 
contour integration methods within the range defined by Ineq. (|148p and, outside this range, by 
means of the method of principal values for singular integrals. 

Fig. 3 Nonequilibrium part of the potential Aipj (+r) is plotted in (x,r) plane at £ = 1. Here 
Aipj (+r) = (+r) / [nze/ \f2irD) with ip® denoting the Debye-Hiickel potential. In Eq. ([92]) ip® 
is not explicitly put in since tpj (±r) is the nonequilibrium part of the potential in the external field. 
Therefore Aipj (+r) should be understood as Aipj (+r) = (+r) — ^ / ^Kze/y27rD^. Within 
the range of x and r satisfying Ineq. (j90|) [also see Ineq. ()A20p ] the contour integration method is 
used and outside the region the method of principal value integration is used for computation. 

Fig. 4 The reduced axial velocity profile v x (x, r, £) is plotted in (x,r) plane at £ = 1. Within 
the range of x and r satisfying Ineq. (I148P the contour integration method is used and outside the 
region the method of principal value integration is used for computation. The axial velocity profile 
is directional, being positive the positive x direction parallel to the external field before vanishing to 
zero at large distance whereas being negative in the transversal (radial) direction before vanishing 
to zero as r increases. Thus the boundary conditions are satisfied in both x and r directions. This 
figure indicates the mode of behaviors of the counterflow of the medium to the ionic movement 
when the external field is turned on. 

Fig. 5 The electrophoretic factor f (x,r, £) is plotted in 3D in a similar color coding to Fig. 4 
in the case of £ = 1. 

Fig. 6 The projection of surface f (x,r, £) onto (x,r) plane. There are two sets of quasi- 
elliptical level curves; one with the major axis on the x axis and the other on the r axis. The 
former corresponds to the contours of the negative part of the f (x, r, £) surface projected onto 
(x, r) plane, and the latter to the contours of the positive part projected onto (x, r) plane. The 
outermost level curve C p is the locus of f (x,r, £) = 0. This level curve Cp depicts the moving 
ion atmosphere distorted by the external electric field from the spherical form assumed by the ion 
atmosphere at £ = 0. This moving ion atmosphere is seen polarized toward the field direction. 

Fig. 7 The distorted ion atmosphere is seen to have the center at (x c , 0) on the x axis. The 
field dependence of the center of the ion atmosphere (x c , 0) describes the trajectory of its motion. 
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The trajectory is shown in this figure. The curve indicates the mode of migration for the center 
from the origin of the coordinate system where the center is located when £ = 0, as the field 
strength is increased. It decreases to a plateau after reaching a maximum as £ increases. 

Fig. 8 Plot of an example for f (x, r, £) at x = r = 0.5 as a function of £ and its comparison 
with Wilson's electrophoretic coefficient /(£)• The solid line, the present theory; the dotted line, 
the OW theory. 

Fig. 9 A 3 — D relaxation time coefficient g (x, r; £). A combination of the contour integration 
results and the method of principal integration is used to construct the surface. 

Fig. 10 Plot of and example for g (x,r,£) at x = r = 0.5 as a function of £ and its comparison 
with Wilson's electrophoretic coefficient g(£). The solid line, the present theory; the dotted line, 
the OW theory. 

Fig. 11 Contour C\ for integrals Kf and Kf . This contour also applies to integrals Jf and Jf 
and P± and P*. The bold line denotes the branch cut. 

Fig. 12 Contour C3 for integrals and K^. This contour also applies to integrals J| and J| 
and P3 and P| . The bold line denotes the branch cut. 

Fig. 13 Contour C2 for integrals K| and Iff. This contour also applies to integrals Jf and Jf 
and Pf an d Pi- The bold line denotes the branch cut on the negative real axis. 
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